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ABSTRACT 


A  theoretical  study  is  presented  for  hydrodynamically  and  therm¬ 
ally  fully  developed  combined  free  and  forced  laminar  convection  with 
upward  flow  in  inclined  tubes  subjected  to  a  uniform  wall  heat  flux. 

The  case  with  uniform  internal  heat  generation  within  the  fluid  is  also 
studied.  A  numerical  solution  using  a  combination  of  boundary  vorticity 
method  and  a  line  iterative  relaxation  method  is  employed,  and  a  com¬ 
parison  is  made  between  the  numerical  analysis  and  the  small  amount  of 
analytical  data  available  in  the  literature.  The  numerical  solution  con¬ 
verges  up  to  a  reasonably  high  value  of  the  characteristic  parameter 
where  generally  an  asymptotic  behavior  for  flow  and  heat  transfer  results 
already  appears,  and  further  results  can  be  obtained  by  a  linear  extra¬ 
polation. 

The  results  of  the  present  investigation  show  that  the  perturbation 
analysis  in  terms  of  power  series  of  Rayleigh  number  has  a  rather  limited 
range  of  applicability  for  the  present  problem,  and  reveal  further  that  a 
maximum  value  for  Mussel t  number  does  not  exist  for  any  tube  inclination 
angle  with  given  values  of  the  dimensionless  parameters  which  is  clearly 
contrary  to  the  result  from  perturbation  solution.  The  tube  inclination 
angle,  or  body-force  orientation  effects  on  flow  and  heat  transfer  charac¬ 
teristics  are  clarified  for  upward  laminar  flow  configuration,  and  show 
that  in  high  Rayleigh  number  regime  the  tube  orientation  effect  has  con¬ 
siderable  influence  on  the  results  in  the  neighborhood  of  horizontal 
direction.  Typical  graphical  results  for  velocity  and  temperature 
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profiles,  flow  and  heat  transfer  parameters  are  presented.  The  influ¬ 
ence  of  each  parameter  such  as  Re,  Ra,  Pr,  F  and  a  on  flow  and  heat 
transfer  results  is  clarified  in  this  study. 
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=  matrix 

=  cross  section  area,  na^  . 

=  radius  of  tube 
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=  axial  temperature  gradient,  3T/3Z 
=  constant  pressure  specific  heat 
=  vector  form 

=  heat  source  parameter,  aQ/p  C  v  C 

P 

”2 

=  friction  factor,  2t  / ( pW  ) ,  or  a  dummy  variable,  see 

w 

equation  (2.12) 

=  vector  form 

=  gravitational  acceleration 
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=  average  heat  transfer  coefficient 
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4 
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CHAPTER  I 

INTRODUCTION 


1 . 1  Review  of  Pertinent  Literature  Relating  to  Laminar  Forced  Convection 

in  Inclined  Tubes  with  Buoyancy  Effect. 

Fully  developed  laminar  convective  heat  transfer  in  various  heated 
or  cooled  ducts  with  double  helix  secondary  flow  caused  by  such  body  for¬ 
ces  as  buoyancy  forces  in  gravitational  or  rotating  field,  centrifugal  or 
Coriolis  forces  acting  in  a  cross-section  normal  to  the  main  flow  has 
been  studied  by  many  investigators  in  recent  years.  It  is  evident  that 
in  many  practical  problems  the  body  forces  cannot  always  be  expected  to 
act  in  a  direction  normal  to  the  main  flow.  It  is  noted  that  the  two 
limiting  cases  of  combined  free  and  forced  laminar  convections  in  verti¬ 
cal  and  horizontal  tubes  have  been  studied  rather  extensively  in  the 
past.  Fully  developed  combined  free  and  forced  laminar  convection  in 
vertical  tubes  or  channels  leads  to  a  linear  partial  differential  equat¬ 
ion  and  various  classical  analytical  methods  can  be  applied  for  the 
solution.  In  contrast  for  the  case  of  combined  free  and  forced  laminar 
convection  in  horizontal  tubes  the  classical  analytical  methods  cannot 
be  readily  applied  since  the  governing  differential  equations  are  non¬ 
linear.  Morton  [1]  approached  the  combined  free  and  forced  laminar 
convection  in  uniformly  heated  horizontal  tubes  at  low  Rayleigh  numbers 
using  a  perturbation  method.  The  perturbation  method  is  of  mathematical 
interest,  but  is  now  known  to  be  applicable  only  for  rather  low  charac¬ 
teristic  parameter  flow  regime  which  is  not  .important  in  design. 
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The  influence  of  tube  orientation  on  fully  developed  combined  free 
and  forced  laminar  convection  was  recently  studied  by  Iqbal  and 
Stachiewicz  [2]  using  a  perturbation  method  identical  to  Morton's  [i]. 

In  view  of  the  fact  that  the  perturbation  method  as  demonstrated  in  the 
literature  for  convective  heat  transfer  with  secondary  flow  diverges 
quickly  with  the  increase  of  the  characteristic  parameter,  it  is  desir¬ 
able  to  study  tube  inclination  angle  or  body-force  orientation  effect  on 
fully  developed  combined  free  and  forced  laminar  convection  in  uniformly 
heated  tubes  by  another  method.  It  is  of  interest  to  point  out  that 
Mori  and  Futagami  [3]  presented  a  theoretical  analysis  for  the  problem 
treated  by  Morton  [1]  which  is  valid  for  high  ReRa  regime  only  by  employ¬ 
ing  boundary-layer  approximation.  In  an  attempt  to  bridge  the  gap  bet¬ 
ween  the  perturbation  method  and  the  boundary-layer  approximation,  Hwang 
and  Cheng  [4]  approached  the  same  problem  numerically  by  using  a  combin¬ 
ation  of  line  iterative  method  and  boundary  vorticity  method.  It  was 
pointed  out  therein  that  the  newly  developed  boundary  vorticity  method  is 
applicable  to  a  general  convective  heat  transfer  with  secondary  flow. 

The  study  of  forced  convection  in  tubes  or  channels  with  internal 
heat  generation  has  also  been  carried  out  by  many  investigators  in  recent 
years  because  of  possible  applications  in  nuclear  and  chemical  reactors. 
Fully  developed  combined  free  and  forced  laminar  convection  in  vertical 
tubes  [5]  or  channels  [6, 7, 8, 9]  with  uniform  internal  heat  generation  has 
been  studied  rather  completely.  In  contrast,  for  the  same  problem  in 
horizontal  tubes,  only  one  study  [10]  using  a  perturbation  method  is  re¬ 
ported  in  the  literature.  In  this  study,  combined  free  and  forced  lami¬ 
nar  convection  in  inclined  tubes  with  uniform  internal  heat  generation 
will  also  be  studied.  The  present  study  can  be  considered  as  an 
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extension  or  generalization  of  the  combined  free  and  forced  laminar  con¬ 
vection  in  vertical  or  horizontal  tubes  reported  in  the  literature. 

1 . 2  Purpose  of  Present  Investigation. 

The  purpose  of  this  investigation  is  to  demonstrate  the  applicability 
of  the  boundary  vorticity  method  [4,11]  for  a  steady  fully  developed  com¬ 
bined  free  and  forced  laminar  convection  in  uniformly  heated  inclined 
tubes  with  and  without  uniform  internal  heat  generation  and  present 
accurate  numerical  results  for  flow  and  heat  transfer.  Besides  clarify¬ 
ing  the  flow  and  heat  transfer  characteristics  in  inclined  tubes, 
comparisons  with  the  available  data  in  the  literature  will  be  made.  In 
order  to  limit  the  scope  of  the  present  study,  consideration  will  be 
given  only  to  the  case  of  fully  developed  upward  laminar  flow  in  the 
inclined  tubes  subjected  to  only  one  type  of  thermal  boundary  condition. 

It  should  be  emphasized  that  the  present  study  represents  a  generalization 
or  extension  of  the  previous  studies  dealing  with  vertical  or  horizontal 
tubes,  and  is  an  asymptotic  case  for  the  Graetz  problem  in  inclined  tubes 
which  apparently  has  not  been  solved  yet. 
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CHAPTER  II 

THEORETICAL  ANALYSIS 


2.1  Formulation  of  the  Problem. 

Consideration  is  given  to  a  steady  hydrodynamically  and  thermally 
fully  developed  forced  laminar  upward  flow  of  viscous  incompressible 
fluid  with  uniform  internal  heat  generation  in  an  inclined  tube  subjected 
to  the  thermal  boundary  conditions  of  axially  uniform  wall  heat  flux  and 
peripherally  uniform  wall  temperature  at  any  axial  position.  The  follow¬ 
ing  assumptions  are  made  to  facilitate  the  theoretical  analysis: 

1.  Physical  properties  of  fluid  are  considered  to  be  constant 
except  density  variation  with  temperature  in  buoyancy  term 
(Boussinesq  approximation). 

2.  Viscous  dissipation  is  negligible. 

3.  Axial  pressure  gradient  is  constant. 

Introducing  cylindrical  coordinates  (R,  <J>,  Z)  as  shown  in  Fig.  1, 
and  applying  the  assumptions  stated  above,  the  governing  equations  for 
the  present  problem  take  the  following  forms: 

Continuity  equation 
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Momentum  equation  in  ^-direction 
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Momentum  equation  in  Z-di recti  on 
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Energy  equation 


p  / 1 ,3T  ,  V  3T  ,  - |3T \  .  /32T 

p  Cp  U3R  R  34)  ^  k^R2 


1  3T  ,  i_ 
R  3R  r2 


32T 

34)2 


)  +  Q 


(2.5) 


Applying  the  no-slip  condition  at  the  wall  and  peripherally  uniform 
wall  temperature  at  any  axial  position,  the  boundary  conditions  for 
equations  (2.2)  to  (2.5)  are: 


U  =  V  =  W  =  T  -  T  =  0  at  pipe  wall 
w  r  r 

It  will  be  assumed  that  the  density  p  in  the  buoyancy  force  term  of  the 
momentum  equation  takes  the  following  form: 

p  “  pw[1+e{Tw  -  T)]’  where  6  = 

In  order  to  express  the  results  in  general  terms, the  following  dimension¬ 
less  dependent  variables,  constants,  parameters,  and  a  dimensionless 
stream  function  are  introduced: 

p 

Dimensionless  radius,  r  =  — 

d 

Dimensionless  velocity, 

U  V  W 

u  "  v/a'  V  '  v/a’  w  Re{v/a) 
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Dimensionless  temperature,  0  = 


VT 

Re  C  a  Pr 


Stream  function. 


Reynolds  number. 


1  U 

11= - X. 

r  34>  5 


v  = 


_ 

8r 


A  a'5  ap 

Re  =  — — «  ,  where  A  =  -  (|y  +  p  g  sin  a) 
4pv  az  w 


Prandtl  number. 


Pr  = 


v 

K 


Rayleigh  number,  Ra  =  -a 


v  K 


Heat  generation  parameter,  F  =  — ^  ^ 


(2.6) 


3T 

where  C  =  ^  =  constant 


The  momentum  and  energy  equations  can  now  be  restated  in  the  following 
dimensionless  forms  after  eliminating  pressure  terms  between  equations 
(2.2)  and  (2.3)  by  cross-differentiation.  The  results  are: 

Momentum  equation  for  secondary  flow 


v2^  =  ufr  +  r  If  +  ReRa(fr  sin  *  +  r  If  cos  ^  cos  “ 


(2.7) 


Vorticity  equation 


Axial  momentum  equation 


=  C 


y2w  =  ufr  +  r  If  +  Ra  0  S1’n  “  '  4 


Energy  equation 


(2.8) 


(2.9) 


=  (“fr  +  r  I?)  ‘  w  +  F 


Re 


(2.10) 


A  -  a  J  q  r  If  '•  •'  ’  -  ■'  ?n  * 


- 


t 
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Through  symmetry  it  is  required  to  consider  only  one-half  of  the 
circular  region  as  shown  in  Fig.  1.  The  boundary  conditions  for 
equations  (2.7)  to  (2.10)  can  be  restated  as  follows: 

at  r  =  1 

(2.11) 

along  <j>  =  0  and  it 

The  governing  equations  (2.7),  (2.9),  (2.10)  are  quasi-! inear  second- 
order  partial  differential  equations  of  elliptic  type,  and  the  vorticity 
equation  (2.8)  is  simply  a  Poisson's  equation.  If  we  substitute  equation 
(2.8)  into  equation  (2.7)  we  obtain  the  governing  equation  for  secondary 
flow  involving  a  biharmonic  operator.  The  resulting  set  of  governing 
equations  for  ip,  w  and  0  was  solved  by  Iqbal  and  Stachiewicz  [2]  using 
a  perturbation  approach  similar  to  Morton's  [1]  for  the  case  without 
heat  generation  in  inclined  tubes  and  by  Iqbal  [10]  with  uniform  internal 
heat  generation  in  horizontal  tubes. 

However,  the  perturbation  method  as  demonstrated  in  the  literature 
for  convective  heat  transfer  with  secondary  flow  is  known  to  diverge 
quickly  with  the  increase  of  the  characteristic  parameter  such  as  ReRa 
for  fully  developed  laminar  flow  in  uniformly  heated  horizontal  tubes. 

In  view  of  the  serious  limitation  of  the  applicability  of  the  pertur¬ 
bation  method,  a  numerical  solution  appears  to  be  the  only  practical 
approach  for  an  accurate  solution  of  the  present  problem.  But  the 
numerical  solution  of  the  biharmonic  equation  in  polar  coordinates  by 
the  available  conventional  method  [12]  using  an  iterative  technique 
with  successive  over  (or  under)  relaxation  is  known  to  converge  extremely 
slowly,  and  is  not  practical  from  the  viewpoint  of  computing  time.  This 


ip  -  =  w  =  0  =  0 


,  _  r  _  30  9w  n 
*  ~  5  '  3?=  0 


* 


8 


difficulty  can  be  overcome  readily  by  using  a  recently  developed  bound¬ 
ary  vorticity  method  [4,11]  for  the  equation  (2.7)  and  (2.8). 


X 

2 . 2  Finite-Difference  Approximation. 

Using  the  three-point  central -difference  operators,  equations  (2.7) 
to  (2.10)  can  be  transformed  into  a  set  of  algebraic  finite-difference 
equations.  The  equations  (2.7)  to  (2.10)  can  be  regarded  as  having  the 
following  form: 


32f  ,  1  3f  L  1  32f  _  ,.3f  ,  V  3f  .  „ 

~2+73F+TrT'u37  +  73?+G 


(2.12) 


3r  '  r“  3(f>‘ 
where  f  is  a  dummy  variable,  u  and  v  are  secondary  velocity  components, 
and  G  represents  the  remaining  terms  in  each  of  the  equations  (2.7),  (2.9) 
and  (2.10).  It  is  noted  that  if  u  =  v  =  0,  equation  (2.12)  reduces  to  a 
Poisson's  equation.  Referring  to  a  coordinate  system  and  numerical  grid 
shown  in  Fig.  1,  the  finite-difference  formulation  of  equation  (2.12) 


can  be  written  as. 


f.  .  .  +  f  _  .  -  2f.  .  f  • « *i  •  —  f .  i  f.  •  i  +  f-  -  2f.  . 
l-l,  J  i+l, J  i,j  +  i+l  ,J  i-l,  J  +  i,j-l  i,J+l  i,J 


(Ar) 


2r.  Ar 


(r.  A<J>)' 


-u  fi+l,j  ~  ,  v  fi ,  j+l  '  fi.j-l  |.  c 

‘  ui,j  2Ar  v1,j  2r .  A<f>  bi,j 


(2.13) 


After  multiplying  by  (Ar)  and  rearranging  equation  (2.13)  becomes 

o  -%  +  r  ui,j)  Vu  -  fi,j  + 

f i+l , j  =  Gi,j  '  ^777^  +  277a^  vi,j]  fi,j-l  '  ^77a^ 


v  if 

2r -A<f>  vi,jJ  Ti,j+1 


(2.14) 


' 
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Letting, 


bi  ■  • 2 

ci  =  1  +  frT  ‘  f"  ui ,  j 


d. 

1 


A 

=  q  -  T  v  If 

i  J  1  ;  Ur.A<J>j  2riA(j>  vi,jJ  fi,j-l 

-  [(tttt)2  -  V,  ,]  f,  4 


r.Acj) 


2r.A<|>  i,jJ  i , j+1 


equation  (2.14)  can  be  reduced  to 


a.  f.  .  •  b  •  "f  •  .  +  c  •  f  • .  *.  .  -  d. 

l  i-l, j  l  i,j  i  i+l, j  l 


(2.15) 


Equation  (2.15)  is  seen  to  be  a  general  finite-difference  expression 
applicable  to  equations  (2.7)  to  (2.10).  In  computing  the  secondary 
velocity  components  using  u  =  ^  anc*  v  =  "  ff-  »  a  more  accurate  result 
can  be  obtained  by  employing  a  five-point  central-difference  formula 
such  as 


3f .  .  , 

i  ,J  


3r 


12Ar  (f5,j  '  6f4, j  +  18f3,j  '  10f2,j  "  for  1  =  2 

TSr  (fi-2,j  '  8fi-l,j  '  8fi+l,j  +  fi+2,j) 

for  i  =  3,4, . . .M  -  1 

W  '  6fM-2,j  +  18fM-l,j  '  10fM,j  '  3fH+l,j)’  for  1  =  M 


'  12A<j>  W,2  +  8fi,3  '  fi,4^’  for  J  2 


(fi,j-2  '  8fi,j-l  +  8fi,j+l  '  fi,j+2)’  for  j  =  3>4>---N-1 


r  J  .  g  **••*>** 
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It  is  noted  that  3f.  ./3<J>  =  0  at  j  =  1  and  N  +  1.  In  order  to  avoid 

"i 

the  singularity  at  the  origin  of  the  cylindrical  coordinates  a  finite- 
difference  equation  in  Cartesian  coordinates  is  employed  at  the  origin. 

2.3  Numerical  Solution  Using  Line  Iterative  Technique  and  Boundary 

Vorticity  Method. 

The  transformation  of  equations  (2.7)  to  (2.10)  into  a  general 
algebraic  finite-difference  equation  (2.13)  was  mentioned  in  Section  2.2. 
The  boundary  conditions  (equation  (2.11))  can  be  written  in  finite- 
difference  forms  as. 


for  j  =  1  and  N  +  1 


(2.16) 


for  i  =  M  +  1 


for  j  =  1  and  N  +  1 


Using  the  above  boundary  conditions,  a  set  of  finite-difference  equations 
may  be  written  in  the  following  forms  after  applying  equation  (2.15)  to 
the  grid  points  along  the  radial  line  j: 


i  =  1 


, . . . 


M-l  (2.17) 


aM  +  bM  fM,j  dM 


i  =  M 


Furthermore,  equation  (2.17)  can  be  written  in  matrix  form  as: 


A  f  =  d 


(2.18) 


' 

11 


or 


bl  C1 

1 

4- 

1 _ 

I 

CL 

_ 1 

a  2  b  2  c  2  0 

f2 

d2 

a3  b3  C3 

•  •  • 

f3 

O 

d3 

• 

•  • 

0  ... 

o 

• 

9 

O 

aM-l  bM-l  CM-1 

fM-l 

dH-l 

aM  bM 

fM 

dM  ‘  CM+1  fM+l 

Here  A  is  seen  to  be  a  tridiagonal  matrix  and  a  Gaussian  elimination 
method  [13]  may  be  used  to  solve  f.  The  following  procedure  using  for¬ 
ward  elimination  and  backward  substitution  is  used, 
h-j  =  c]/b1 

h .  =  c./(b.  -  a.h.  -j)  i  =  2,3, . M-l 

/ 

q1  =  d1/b1 

q.  =  (d.  -  a.  q^p  /  (tK  -  a.h.^)  1  =  2,3 . M 

and  fM  = 

fi  =  qi  *  hifi+l  1  =  M'1, 


It  is  noted  that  the  numerical  solution  can  be  applied  to  equat¬ 
ions  (2.9)  and  (2.10)  with  the  boundary  conditions  given  in  equation 
(2.16),  but  the  boundary  condition  for  equation  (2.7)  is  unknown.  On 
the  other  hand,  two  boundary  conditions  are  available  for  equation  (2.8). 
The  success  of  the  boundary  vorticity  method  is  based  on  the  observation 
that  a  linear  relationship  exists  between  the  vorticity  £M+]  .  and  the 


' 


/■')  *  jl> 
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stream  function  ^M+1  .  at  the  boundary  [4.11].  Using  the  boundary 
condition  of  i|k_i  j  ^  j+-j  j  for  i  =  M+l ,  some  elements  of  the  matrix 
equation  (2.18)  such  as  a^,  b^,  f^  and  ^"cjv]+i^[vi+-]  W1*^  now  be  replaced 

by  aM+i+cM+i *  bM+ls  fM+l  and  dM+ls  respectively.  Assuming  arbitrary 
values  and  for  ^+-j  .  and  solving  finite-difference  equations 
for  equations  (2.7)  and  (2.8),  the  values  for  4^  and  ip^  can  be 
obtained  using  the  boundary  condition  |^  =  0  at  wall .  The  boundary 
values  for  the  vorticity  £M+1  j  and  the  stream  function  .  are 
related  by  the  following  linear  equation: 


M)  J 1) 

■(3)  -  h>  ~  ,  (3)  ,„(2) 

’b  7(2i  rnr 


+  5 


(2) 


(2.19) 


Using  this  linear  relation  and  noting  that  the  correct  stream  function 

(3)  (3) 

at  the  boundary  is  ^5  =  0,  ££  can  be  computed.  With  the  boundary 

vorticity  numerically  determined,  the  numerical  solution  of  a  set  of 
equations  (2.7)  to  (2.10)  can  be  carried  out  by  a  line  iterative  relax¬ 
ation  procedure  using  the  following  equation: 


f  (n+l ) 


+  “  [fi,j 


(n+l) 


(2.20) 


where  w  is  a  relaxation  factor.  The  numerical  solution  of  a  set  of 

equations  (2.7)  to  (2.10)  can  be  started  by  assuming  initial  values  for 

w.  6.  u.  .,  and  v.  .  with  the  initial  mesh  size  of  M,N  =  14.  The 
i,j  i,j  i  ,j  i  ,J 

numerical  solution  will  be  continued  until  the  dependent  variables 

w.  .,6.  .,  ij;.  .,  and  £.  .  satisfy  the  following  prescribed  error: 
i  jj  i  >J  i  >J  i  »J 


The  secondary  velocity  components  u.  .  and  v.  .  can  be  obtained  by  using 

1  jJ  ■ 


' 
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the  five-point  finite-difference  formula  for  the  first  derivatives. 

In  order  to  obtain  more  accurate  numerical  results,  the  procedure  is 
repeated  by  using  the  mesh  size  of  M,N  =  28  and  iteration  is  continued 
until  the  prescribed  error  is  reached. 

The  method  of  determining  an  optimum  relaxation  factor  is  not 
available  for  the  present  non-linear  problem.  However,  it  is  found  that 
the  influence  of  overrelaxation. (w  >  1)  or  underrelaxation  (w  <1)  on 
computing  time  is  not  appreciable.  Consequently,  the  relaxation  factor 
of  unity  is  used  in  most  of  the  calculations.  The  rate  of  convergence 
for  line  iterative  method  depends  on  the  values  of  the  parameters  Re, 

Ra,  Pr,  a  and  F.  For  given  values  of  the  parameters  Re,  Ra,  Pr,  a  and 
F,  it  is  possible  to  obtain  flow  and  heat  transfer  results  within  3 
minutes  in  the  low  Rayleigh  number  regime  on  the  IBM360/67  system.  How¬ 
ever,  in  high  parameter  regime  such  as  Ra  >  4,000,  it  takes  around  5 
minutes  to  complete  one  calculation  for  fixed  values  of  the  parameters 
Re,  Ra,  Pr,  a  and  F.  It  is  noted  that  in  very  high  Rayleigh  number 
regime  the  matrix  cannot  be  inverted.  In  order  to  check  the  accuracy  of 
the  numerical  results  two  different  methods  of  obtaining  the  Nusselt 
number  and  the  pressure-drop  parameter  are  used.  Since  the  prescribed 
error  e  is  set  to  be  less  than  10“5,  it  is  expected  that  the  numerical 
results  for  the  pressure  drop  parameter  and  Nusselt  number  will  be  valid 
up  to  four  significant  figures.  This  observation  is  confirmed  by  the 
agreement  between  the  numerical  results  based  on  two  alternative  defin¬ 
itions  for  the  flow  and  heat  transfer  results.  Typical  numerical 
results  are  given  in  Appendix  C. 


- 


14 


CHAPTER  III 

FLOW  AND  HEAT  TRANSFER  RESULTS  WITHOUT  HEAT  SOURCES 
3.1  Velocity  Distributions. 

For  the  present  problem  without  heat  sources,  four  independent 
parameters  Re,  Ra,  Pr,  and  a  appear  in  the  governing  equations,  and  it 
is  not  practical  to  carry  out  numerical  computation  for  all  possible 
cases.  Consequently,  consideration  will  be  limited  to  a  few  typical 
cases  to  illustrate  the  effect  of  each  parameter. 

The  axial  velocity  profile  depends  on  Prandtl  number,  Rayleigh  num¬ 
ber,  Reynolds  number,  and  inclination  angle.  The  effect  of  each  para¬ 
meter  on  axial  velocity  profile  along  <j>  =  0  and  tt  is  illustrated  in 
Fig.  2(a)  to  (d).  The  effect  of  Rayleigh  number  for  Pr  =  10,  Re  =  20, 
and  a  =  45°  demonstrated  in  Fig.  2(a)  shows  that  at  Ra  =  1,000,  back 
flow  already  occurs  from  r  =  0  to  r  =  0.25  along  <|>  =  it.  The  Rayleigh  num¬ 
ber  at  which  the  back  flow  occurs  is  of  considerable  theoretical  and 
practical  interest,  and  the  numerical  result  is  given  below  for  Pr  =  0.75 
and  ReRa  _  50,000. 

a  0°  10°  20°  30°  45°  60°  75°  90° 

Ra  oo  3,800  1,950  1,350  930  780  690  670 

It  is  noted  that  for  pure  forced  convection  with  Ra  =  0  the  axial 
velocity  profile  is  parabolic.  As  the  Rayleigh  number  increases,  the 
intensity  of  the  secondary  motion  increases.  This  secondary  motion  is 
known  to  distort  the  parabolic  velocity  profile,  and  decreases  the  mag¬ 
nitude  of  the  maximum  velocity. 

The  effect  of  Reynolds  number  appears  through  a  parameter  ReRa  in 
the  momentum  equation  (2.7)  for  secondary  flow,  and  the  role  of  Reynolds 


15 


number  on  axial  velocity  profile  is  seen  to  be  similar  to  that  of  ReRa 
for  the  limiting  case  of  horizontal  tube.  Fig.  2(b)  shows  that  the  maxi¬ 
mum  axial  velocity  moves  toward  the  direction  cf>  =  tt  as  the  Reynolds 
number  increases. 

For  given  values  of  Pr,  Re  and  Ra,  the  effect  of  inclination  angle 
a  can  be  seen  from  axial  momentum  equation  (2.9)  using  membrane  analogy, 
and  Fig.  2(c)  confirms  the  well  known  axial  velocity  profiles  for  the 
horizontal  (a  =  0°)  and  vertical  tubes  (a  =  90°). 

The  effect  of  Prandtl  number  appears  in  the  convective  terms  of  the 
energy  equation  (2.10).  As  the  Prandtl  number  increases,  the  secondary 
motion  decreases,  and  the  axial  velocity  profile  approaches  that  of  Hagen- 
Poiseuille  flow,  as  shown  in  Fig.  2(d).  This  observation  is  significant 
since  it  leads  to  the  possibility  for  the  numerical  solution  of  Graetz 
problem  with  significant  buoyancy  effect.  Fig.  2(d)  shows  the  effect  of 
Prandtl  number  on  axial  velocity  profile  for  given  values  of  ReRa,  Ra 
and  a.  The  magnitude  of  the  maximum  velocity  is  seen  to  increase  as 
Prandtl  number  increases. 

Typical  streamlines  are  illustrated  in  Fig.  3  for  a  =  45°  and  the 
patterns  are  seen  to  be  similar  to  those  for  horizontal  tube  a  =  0°  [4]. 
There  are  two  symmetric  streamline,  patterns  for  the  secondary  flow.  One 
notes  that  at  low  Rayleigh  number  the  centers  of  circulation  first 
appear  on  the  radial  line  <j>  =  or  -  With  the  increase  of  the 
Rayleigh  number,  they  tend  to  move  toward  the  wall  radially,  but  circum¬ 
ferentially  they  move  away  from  (j>  =  or  -  j  toward  <j>  =  tt.  With  further 
increase  of  Rayleigh  number,  they  then  move  back  toward  the  radial  line 
<j>  =  ~  or  -  For  given  values  of  Pr,  Re  and  Ra,  the  location  of  the 
eye  of  the  streamlines  does  not  appear  to  be  sensitive  to  the  inclination 
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angle  effect,  but  the  absolute  value  of  the  stream  function  at  the  eye 
decreases  with  the  increase  of  the  inclination  angle.  Of  course,  the 
secondary  flow  disappears  completely  at  a  =  90°.  The  distributions  of 
the  secondary  flow  velocity  components  u,  v,  will  be  discussed  later 
in  Chapter  IV. 

3.2  Temperature  Distributions. 

The  effects  of  the  parameters  Ra,  Re,  a,  and  Pr,  respectively,  on 
the  temperature  profile  along  <J>  =  0  and  7r  are  shown  in  Fig.  4(a)  to  (d) 
under  the  corresponding  conditions  as  those  for  the  axial  velocity 
profiles  discussed  in  Section  3.1.  Qualitatively,  the  profiles  for  tem¬ 
perature  and  axial  velocity  are  similar  except  the  phenomenon  of  reverse 
flow.  This  observation  is  confirmed  by  the  similarity  of  the  respective 
governing  equation. 

Fig.  4(a)  shows  that  at  Ra  =  1,000  with  back  flow  (see  Fig.  2(a)), 
the  secondary  flow  is  predominant  and  the  temperature  profile  is  flat  in 
the  central  core  region,  indicating  that  the  conduction  terms  are  neglig¬ 
ible  as  compared  with  the  advective  terms  in  the  energy  equation  (2.10). 
Fig.  4(b)  shows  that  the  effect  of  Reynolds  number  on  temperature  profile 
is  qualitatively  similar  to  that  of  Reynolds  number  effect  on  axial  velo¬ 
city  profile.  As  Reynolds  number  increases,  the  magnitude  of  the  maximum 
temperature  difference  6  decreases,  and  its  location  shifts  toward  the 
lower  wall.  In  addition,  the  temperature  gradient  at  the  bottom  is  seen 
to  be  higher  than  that  at  the  top  of  the  tube  for  high  Reynolds  number 
regime.  For  given  values  of  the  parameters  Re,  Ra,  and  Pr,  the  effect  of 
inclination  angle  on  temperature  profile  is  shown  in  Fig.  4(c).  For  the 
limiting  case  of  vertical  tube,  the  temperature  profile  is  symmetric.  In 
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other  words,  the  isotherms  are  concentric  for  the  vertical  tube  case. 

The  concentric  isotherms  will  be  gradually  distorted  as  the  tube  inclin¬ 
ation  angle  a  varies  from  a  =  90°  to  a  =  0°,  and  the  center  of  the  iso¬ 
therms  will  move  downward  along  the  direction  <J>  =  tt.  For  given  values  of 
Re,  Ra  and  a,  the  effect  of  Prandtl  number  on  the  temperature  profile  is 
qualitatively  similar  to  that  of  Ra  or  Re  as  can  be  seen  from  the  energy 
equation  (2.10).  This  observation  is  also  confirmed  in  Fig.  4(d).  One 
sees  that  as  Prandtl  number  increases,  the  symmetric  temperature  distri¬ 
bution  is  distorted  gradually,  and  the  magnitude  of  the  maximum  temper¬ 
ature  difference  decreases.  Comparing  the  effects  of  Prandtl  number  on 
axial  velocity  and  temperature  profiles,  it  is  clear  that  the  effects  of 
Prandtl  number  on  axial  velocity  and  temperature  profiles  are  seen  to  be 
completely  opposite. 

The  typical  isotherms  are  illustrated  in  Fig.  3  for  a  =  45°,  and  the 
pattern  is  seen  to  be  similar  to  that  for  horizontal  tube  [4].  The  tem¬ 
perature  drops  rather  gradually  from  the  value  at  the  wall  to  the  minimum 
value  along  the  direction  <f>  =  0.  Along  the  direction  <j>  =  tt,  the  temper¬ 
ature  drops  rapidly  from  the  wall  temperature  to  the  same  minimum  value. 

A  study  of  the  isotherms  in  high  Rayleigh  number  regime,  for  example  at 
Ra  =  4,000  with  a  =  45°,  reveals  that  the  isotherms  are  nearly  circular 
and  concentric.  This  situation  can  also  be  seen  from  the  temperature 
profile  for  Ra  =  1,000  shown  in  Fig.  4(a),  and  is  apparently  related  to 
back  flow  phenomenon  shown  in  Fig.  2(a).  Noting  that  the  isotherms  are 
concentric  for  the  case  without  secondary  flow,  one  sees  that  as  the 
Rayleigh  number  increases  from  zero,  the  concentric  isotherms  will  be  grad¬ 
ually  distorted  and  finally  become  concentric  again  at  sufficiently  large 
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Rayleigh  number.  It  is  noted  that  the  pattern  of  the  isotherms  is  close¬ 
ly  related  to  the  movement  of  the  centers  of  circulation  for  the  stream¬ 
lines.  When  the  centers  of  circulation  return  to  the  radial  line  <j>  =  -| 
or  -  y  again  at  sufficiently  large  Rayleigh  number,  the  isotherms  become 
nearly  concentric  circles.  One  should  note  that  in  discussing  the 
Rayleigh  number  effect,  the  remaining  parameters  such  as  Pr,  a  and  Re  are 
held  constant. 


3.3  The  Effects  of  Ra  and  Re  on  Flow  and  Heat  Transfer  Results. 


The  expressions  for  the  product  of  friction  factor  and  Reynolds  num¬ 
ber,  f*Re,  and  Nusselt  number,  Nu,  can  be  obtained  by  considering  either 
the  velocity  and  temperature  gradients,  respectively ,  along  the  tube  wall 
or  the  overall  force  and  energy  balances,  respectively ,  for  the  axial 
length  dZ. 

Following  the  usual  definition  for  the  product  of  friction  factor 
and  Reynolds  number,  we  obtain 


f-Re  = 


Tw  .  W  2a  p  =  4a  Tw 


W2 


u 


Ji  w 


(3.1) 


where  tw  can  be  obtained  by  considering  the  velocity  gradient  along  the 


tube  wall  as. 


Tw 


.  (f£)  R  d*  _ 

s  3R  W  _  .  (3W) 


R  d<J> 


w 


substituting  equation  (3.2)  into  equation  (3.1),  we  have 


(f  Re ) -j  = 


4  (— ) 


w 


w 


(3.2) 


(3.3) 
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A  force  balance  on  a  fluid  element  Ar  dZ  gives 

\  '  S  =  -  Ar  (ll  +  P  9  s1n  «) 

=  Ar  A ( 1  -  1  Ra  9  sin  a)  (3.4) 

Substituting  equation  (3.4)  into  equation  (3.1)  and  using  the  dimension¬ 
less  mean  velocity  w,  we  obtain 


(f*Re)?  =  ~  [  1  -  Ra  6  ]  (3.5) 

L  w  ^ 

In  the  limiting  case  of  horizontal  tube  a  =  0,  equation  (3.5)  reduces  to 


(t  Re)?  =  I 
L  w 

Analogous  to  the  resistance  coefficient,  the  Nusselt  number  can  be 
obtained  in  two  ways.  One  way  is  to  consider  the  temperature  gradient  at 
the  wall,  and  the  other  way  is  to  consider  the  overall  energy  balance  for 
the  axial  length  dZ.  The  Nusselt  number  is  defined  as 

Nu  =  (3.6) 


Using  the  temperature  gradient  at  wall,  the  average  heat  transfer 
coefficient  h  can  be  written  as. 


or 


where 


h  TM  S  dZ 


=  -k  (|I)  S  dZ 


W 


h  =  - 


w 


M 


tm  = 


W  (Tw  -  T)  d  Ar 


U  d  A. 


=  Re  C  a  Pr 


W0 


w 
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So,  the  Nusselt  number  can  be  expressed  as, 


2-w(f) 

w 


w 


(3.7) 


W0 


On  the  other  hand,  by  considering  the  overall  energy  balance,  we  have 


pCpl“Ar^THS 


or 


p  Cp  C  Ar  Re  v 


r 


h  = 


w 


Consequently,  equation  (3.6)  can  be  expressed  as 

_2 

(Nu)?  =  ^  (3.8) 

we 

For  the  present  problem  without  heat  sources,  the  parameters  ReRa, 

Ra,  Pr  and  a  appear  in  the  governing  equations,  and  each  parameter  has  its 
influence  on  flow  and  heat  transfer  results.  Fig.  5  shows  the  ratio 
fRe/(fRe)g  as  a  function  of  Rayleigh  number  with  inclination  angle  a  as  a 
parameter  for  Pr  =  0.75  and  ReRa  =  4,000.  It  is  seen  that  after  reaching 
approximately  a  value  of  2  for  fRe/(fRe)Q,  an  asymptotic  behavior  appears 
for  each  curve,  and  all  the  curves  are  seen  to  be  parallel  to  each  other. 
As  the  Rayleigh  number  increases,  secondary  motion  increases,  so  does  the 
fRe/(fRe)Q.  Particularly  at  high  Rayleigh  number  regime,  the  product  fRe 
increases  considerably  as  a  varies  from  0°  to  45°.  One  notes  that  a 
slight  deviation  from  the  horizontal  direction  a  =  0°  introduces  a  sub¬ 
stantial  increase  in  friction  factor  in  the  high  Rayleigh  number  regime. 
With  Pr  and  a  given,  it  is  seen  that  the  value  of  fRe  depends  on  Ra  only 
after  reaching  a  certain  value  of  Ra.  This  effect  is  clearly  seen  in 
Fig.  6,  where  the  variation  of  fRe/(fRe)Q  with  Ra  is  shown  for  Pr  =  0.75 
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and  a  =  45°  with  ReRa  as  a  parameter.  After  reaching  say  Ra  =  200,  the 
value  for  fRe  is  seen  to  be  independent  of  Re  and  depends  only  on  Ra. 

Note  that  the  limiting  value  for  Re  is  about  2,000  for  laminar  flow. 

The  effects  of  Rayleigh  number  on  heat  transfer  results  and  flow 
results  are  known  to  be  similar,  and  the  heat  transfer  results  are  shown 
in  Fig.  7  and  Fig.  8.  Fig.  7  shows  that  after  reaching  a  certain  value 
of  Ra,  the  Nusselt  number  ratio  Nu/(Nu)q  is  a  function  of  Ra  only  for 
each  a.  The  asymptotic  behavior  in  the  form  of  parallel  straight  lines 
for  various  a's  in  the  high  Rayleigh  number  regime  as  shown  in  Fig.  7 
suggests  that  a  new  correlation  is  possible  for  heat  transfer  results,  and 
this  is  demonstrated  in  Fig.  8  where  all  the  heat  transfer  results  in 
Fig.  7  are  now  replotted  on  the  basis  of  Nu/(Nu)Q  versus  Ra  sin  a.  After 
reaching  a  value  of  say  Ra  sin  a  =  70,  all  the  heat  transfer  results  for 
the  range  0  <  a  1  90°  can  be  predicted  by  a  single  curve  with  sufficient 
accuracy.  The  tube  orientation  effect  for  Ra  sin  a  <  70  is  of  interest. 

The  effect  of  Reynolds  number  on  heat  transfer  results  is  shown  in 
Fig.  9  to  Fig.  12.  Fig.  9  shows  that  after  reaching  Ra  =  450,  the 
Nusselt  number  ratio  Nu/(Nu)Q  is  independent  of  Reynolds  number,  and  is  a 
function  of  Ra  only.  The  curves  for  ReRa  =  50,000  and  10,000  show  that 
minimum  value  for  Nu/(Nu)Q  exists  at  some  value  of  Ra.  An  examination  of 
the  momentum  equation  (2.7)  and  the  energy  equation  (2.10)  shows  that  for 
horizontal  tube  case  (a  =  0°),  the  flow  and  heat  transfer  results  depend 
on  Pr  and  ReRa.  On  the  other  hand,  for  the  vertical  case  (a  =  90°),  the 
flow  and  heat  transfer  results  are  independent  of  Re  and  Pr.  This  suggests 
that  the  Reynolds  number  effect  is  significant  in  the  neighborhood  of 
horizontal  case.  Figs.  10  to  12  show  the  effect  of  Reynolds  number  on 
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Nu/(Nu)q  for  a  =  0°,  30°,  and  60°»  respectively ,  with  Pr  =  5.  It  is  clear 
that  for  the  horizontal  case  shown  in  Fig.  10,  the  heat  transfer  rate 
increases  with  the  increase  of  Reynolds  number  for  a  given  value  of  Ra. 
Furthermore,  for  various  values  of  Re,  the  curves  are  seen  to  be  parallel 
to  each  other  in  the  high  Rayleigh  number  regime.  Two  curves  for  Re  =  5 
and  Re  =  20  in  the  case  of  a  =  30°  and  Pr  =  5  shown  in  Fig.  11  start  sig¬ 
nificant  deviation  at  around  Ra  =  10.  This  indicates  that  the  Reynolds 

number  effect  appears  after  reaching  around  Ra  =  10.  However,  after 
reaching  say  Ra  =  1,000,  the  effect  of  Reynolds  number  is  seen  to  be  neg¬ 
ligible.  Fig.  12  shows  that  at  a  =  60°  the  influence  of  Reynolds  number 

on  heat  transfer  results  is  not  appreciable  for  any  Rayleigh  number.  It 
is  seen  that  two  curves  for  Re  =  5  and  Re  =  20  in  Fig.  12  coincide  after 
reaching  around  Ra  =  700.  A  comparison  between  Figs.  10  and  12  shows 
clearly  that  the  effect  of  Reynolds  number  on  heat  transfer  result  is  sig¬ 
nificant  for  the  horizontal  tube  (a  =  0°)  but  insignificant  at  a  =  60°. 

It  is  expected  that  from  a  =  60°  up  to  the  limiting  vertical  case  (a  =  90°) 
the  heat  transfer  result  will  be  independent  of  Reynolds  number. 

The  numerical  results  from  this  study  shown  in  Figs.  10  to  12  can  be 
compared  against  the  results  from  perturbation  solution  [14].  Figs.  10 
to  12  show  clearly  that  there  is  a  significant  difference  between  the  pre¬ 
sent  numerical  solution  and  the  perturbation  solution.  It  should  be 
pointed  out  that  for  the  horizontal  tube  case  the  perturbation  solution 
is  known  to  diverge  quickly  with  the  increase  of  the  parameter,  and  the 
rather  sharp  increase  in  the  Nusselt  number  with  the  increase  of  Ra  shown 
in  Fig.  10  for  cx  =  0°  is  typical  of  the  perturbation  solution.  It  is 
noted  that  for  horizontal  tube,  the  Nusselt  number  is  dependent  on  the 
parameter  ReRa.  The  trend  of  heat  transfer  results  from  perturbation 
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analysis  for  the  inclination  angles  a  =  30°  and  60°  shown  in  Fig.  11  and 
Fig.  12  respectively,  also  shows  the  same  blow-up  trend  similar  to  that 
for  a  =  0°. 

3.4  The  Effects  of  Prandtl  Number  and  Tube  Inclination  Angle  on  Flow 

and  Heat  Transfer  Results. 

The  Prandtl  number  effects  on  flow  and  heat  transfer  results  for 
a  =  45°  and  ReRa  =  4,000  are  shown  in  Fig.  13  and  Fig.  14  respectively. 

Fig.  13  shows  that  the  Prandtl  number  effect  on  fRe/(fRe)g  is  signific¬ 
ant  only  in  the  low  Rayleigh  number  regime.  As  Prandtl  number  increases, 
the  ratio  fRe/(fRe)g  is  seen  to  decrease  in  the  low  Rayleigh  number  region. 
After  reaching  about  Ra  =  2,000,  the  ratio  fRe/(fRe)Q  is  seen  to  be  indep¬ 
endent  of  Prandtl  number.  Fig.  14  shows  Nu/(Nu)q  versus  Ra  with  Pr  as  a 
parameter.  The  Prandtl  number  effect  is  significant  in  the  low  Rayleigh 
number  regime.  Note  the  completely  opposite  effects  of  Prandtl  number  on 
flow  and  heat  transfer  results  in  the  low  Rayleigh  number  regime.  It  is 
noted  that  Prandtl  number  appears  in  the  convective  terms  of  the  energy 
equation  (2.10).  Qualitatively,  the  effect  of  Prandtl  number  on  heat 
transfer  rate  is  similar  to  that  of  secondary  motion.  In  other  words, 
heat  transfer  rate  increases  with  the  increase  of  the  Prandtl  number  in 
the  low  Rayleigh  number  region. 

The  effect  of  tube  inclination  angle  on  flow  and  heat  transfer  res¬ 
ults  is  of  considerable  interest,  and  is  shown  in  Figs.  15  to  18.  It 
is  seen  that  in  the  high  Rayleigh  number  regime  a  slight  variation  of 
inclination  angle  introduces  a  substantial  increase  in  fRe/(fRe)g  and 
Nu/(Nu)q.  The  practical  implication  of  this  effect  is  that  the  inclin¬ 
ation  angle  effect  cannot  be  neglected,  for  example,  for  helical  tubes 
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in  high  Rayleigh  number  regime. 

The  variation  of  fRe/(fRe)Q  with  a  is  shown  in  Fig.  15  for  Pr  =  0.75 
and  ReRa  =  10,000  with  Ra  as  a  parameter.  At  Ra  =  10,  tube  orientation 
for  the  range  a  =  0°  ~  30°  is  seen  to  be  insensitive  to  the  value  of  fRe, 
while  at  Ra  =  4,000,  tube  inclination  effect  on  the  value  of  fRe  is  con¬ 
siderable  in  the  range  of  a  =  0°  -20°.  In  Fig.  16,  the  ratio  Nu/(Nu)Q 
is  plotted  against  a  with  Ra  as  a  parameter  for  Pr  =  0.75  and  ReRa  = 

10,000.  It  is  seen  that  the  inclination  angle  effect  on  Mussel t  number 
ratio  in  low  Rayleigh  number  regime  is  opposite  to  that  in  the  high 
Rayleigh  number  regime.  The  curves  for  Ra  =  400  and  100  show  that  mini¬ 
mum  value  for  Nu/(Nu)q  exists  at  some  value  of  a.  At  high  Rayleigh  num¬ 
ber,  say  Ra  =  4,000,  heat  transfer  result  changes  significantly  with 
slight  deviation  from  horizontal  direction  a  =  0°  which  suggests  again 
that  inclination  angle  effect  cannot  be  neglected  for  certain  helical 
tubes. 

The  problem  under  consideration  was  solved  by  Iqbal  and  Stachiewicz 
[2]  using  perturbation  method  and  it  is  desirable  to  compare  the  results 
from  this  work  with  those  of  references  [2]  and  [14].  Fig.  17  shows  that 
there  is  a  substantial  difference  between  the  present  numerical  solution 
and  the  analytical  solution  using  perturbation  method.  The  present  numer¬ 
ical  solution  checks  exactly  with  the  known  numerical  results  for  hori¬ 
zontal  tube  (a  =  0°)  [4]  and  the  vertical  tube  (a  =  90°)  [5],  whereas  the 
perturbation  solution  [2]  covers  the  range  10°  <  a  -  50°  only.  The  curves 
for  ReRa  =  3,000,  2,000,  and  1,000  from  perturbation  solution  do  not 
appear  to  approach  the  value  for  a  =  90°.  In  view  of  the  considerable  dis¬ 
crepancy  noted  above,  it  is  believed  that  the  perturbation  solution  shown 
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in  reference  [2]  is  invalid.  In  particular,  the  suggestion  [2]  that  for 
any  combination  of  Rayleigh,  Reynolds,  and  Prandtl  numbers,  there  is  an 
optimum  value  of  tube  inclination  which  gives  the  maximum  value  of 
Nusselt  number  for  a  ranging  from  20  to  60  degrees  appears  to  be  erron¬ 
eous.  Fig.  18  further  demonstrates  the  divergence  of  the  perturbation 
solution.  Once  again,  the  suggestion  of  reference  [14]  that  there  is  a 
maximum  heat  transfer  result  in  the  range  of  20  to  60  degrees  of  inclin¬ 
ation  angle  appears  to  be  incorrect.  Fig.  18  shows  that  even  at  Re  = 

4,  6,  8,  the  perturbation  solution  diverges  for  Pr  =  5  and  Ra  =  100.  It 
is  reasonable  to  conclude  that  the  perturbation  solution  for  the  present 
problem  is  valid  only  for  very  low  Rayleigh  number  regime,  say  up  to 
Ra  5  20,  and  there  is  no  optimum  inclination  angle  for  any  given  values 
of  the  Rayleigh,  Reynolds,  and  Prandtl  numbers. 
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CHAPTER  IV 

FLOW  AND  HEAT  TRANSFER  RESULTS  WITH  UNIFORM  HEAT  SOURCES 

4 . 1  Influence  of  Heat  Sources  on  Velocity  and  Temperature  Profiles. 

The  axial  velocity  and  temperature  profiles  depend  on  Prandtl 
number,  tube  inclination  angle,  Reynolds  number,  Rayleigh  number  and 
heat  generation  function.  Since  the  effect  of  each  parameter  on  axial 
velocity  and  temperature  distributions  is  illustrated  in  Chapter  III  for 
the  case  without  heat  generation,  for  brevity,  only  the  effects  of  heat 
generation  and  Rayleigh  number  will  be  considered  here. 

The  effect  of  dimensionless  heat  generation  parameter  F  on  the  dimen¬ 
sionless  axial  velocity  and  temperature  profiles  along  <j>  =  0  and  tt  is 
shown  in  Fig.  19(a)  and  (b),  respectively,  for  Pr  =  0.75,  a  =  45°,  Re  = 
100,  and  Ra  =  40.  It  can  be  seen  that  the  velocity  distribution  near 
the  wall  becomes  steeper  with  larger  values  of  F.  This  fact  agrees,  for 
example,  with  the  result  for  combined  free  and  forced  convection  with 
uniform  internal  heat  generation  along  a  vertical  flat  plate  reported  in 
the  literature.  The  shift  of  the  location  of  the  maximum  velocity  toward 
the  direction  <j>  =  ir  shown  in  Fig.  19(a)  for  a  given  value  of  F  is  known 
to  be  caused  by  secondary  flow.  As  the  value  of  the  heat  sources  para¬ 
meter  F  increases,  the  maximum  velocity  increases  and  the  location  of  the 
maximum  velocity  tends  to  shift  back  toward  the  pipe  axis.  It  is  expected 
that  at  a  sufficiently  high  value  of  F,  the  velocity  distribution  will  be 
parabolic.  This  trend  is  expected  since  with  the  increase  of  heat  sour¬ 
ces  the  temperature  difference  between  the  wall  and  fluid  will  decrease 
(see  Fig.  19(b)).  Noting  the  above  observations,  the  increase  of  the 
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dimensionless  axial  velocity  with  increasing  value  of  F  can  be  seen 
qualitatively  from  the  axial  momentum  equation  (2.9).  Fig.  19(b)  shows 
that  the  dimensionless  temperature  difference  0  is  lower  near  the  wall 
and  elsewhere  with  larger  values  of  F.  The  shift  of  the  location  of  the 
maximum  temperature  difference  is  similar  to  that  for  the  velocity  pro¬ 
file  and  can  be  explained  from  the  free  convection  or  secondary  flow 
effect.  At  F  =  50,  it  is  seen  that  the  temperature  difference  0  becomes 
negative  all  around  the  wall  (r  =  1)  indicating  the  reversed  heat  flow 
there.  The  behavior  of  the  temperature  profiles  shown  in  Fig.  19(b)  can 
also  be  seen  from  the  relative  contribution  of  the  heat  sources  term  and 
the  axial  convection  term  on  the  right  hand  side  of  the  energy  equation 

(2.10),  keeping  in  mind  that  the  secondary  flow  decreases  with  the 

« 

increase  of  heat  sources  parameter  F. 

The  effect  of  Rayleigh  number  on  the  axial  velocity  and  temperature 
profiles  along  <f>  =  0  and  tt  is  shown  in  Fig.  20(a)  and  (b)  respectively, 
and  is  seen  to  be  similar  to  the  case  without  heat  sources  F  =  0  which 
is  shown  in  Chapter  III.  The  trend  of  the  Rayleigh  number  effect  on  the 
axial  velocity  profile  shown  in  Fig.  20(a)  suggests  that  at  sufficiently 
high  Rayleigh  number  the  back  flow  may  occur.  At  Ra  =  1,000  the  temper¬ 
ature  distribution  0  is  already  negative  near  the  wall,  indicating  rever¬ 
sed  heat  flow  situation.  Of  course  this  phenomenon  cannot  occur  without 
heat  sources.  In  Fig.  20(b)  the  shift  of  location  and  the  magnitude  of 
the  maximum  0  is  apparently  associated  with  the  intensity  of  secondary 
flow.  Fig.  20(a)  also  shows  that  the  distortion  from  parabolic  profile 
increases  with  the  increase  of  Rayleigh  number.  The  effects  of  heat 
sources  parameter  and  Rayleigh  number,  respectively  on  the  velocity  and 
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temperature  profiles  shown  in  Fig.  19  and  Fig.  20  with  all  other  para¬ 
meters  being  held  constant  will  be  seen  to  be  useful  in  interpreting 
the  flow  and  heat  transfer  results. 

In  view  of  the  fact  that  the  secondary  flow  velocities  are  too  small 
to  be  measured  in  experiment,  the  distributions  of  the  secondary  velocity 
components  u,  v  are  of  some  interest,  and  these  are  shown  in  Fig.  21(a) 
and  (b)  for  several  representative  values  of  Ra  with  a  -  45°,  Pr  =  0.75, 
Re  =  100  and  F  =  20.  The  trend  of  the  effect  of  Rayleigh  number  on  sec¬ 
ondary  velocity  distribution  indicates  that  boundary  layer  approximation 
does  not  appear  to  be  possible  even  at  high  Rayleigh  number  for  F  =  20. 


4 . 2  Heat  Sources  Effect  on  Flow  and  Heat  Transfer  Results. 

The  expressions  for  the  product  of  friction  factor  and  Reynolds  num¬ 
ber  (fRe)  and  the  Nusselt  (Nu)  are  the  same  as  those  given  in  Chapter  III 
for  the  case  without  heat  generation  except  that  the  second  definition 
for  the  Nusselt  number  must  be  modified.  From  the  energy  balance  for  the 
fluid  element  Ar  dZ,  we  obtain 
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P  Cp  C  Wc  w  Ar  -  Q  Ar  =  h  T^  S, 
where  W  =  Re  (v/a) 

V 

The  average  heat  transfer  coefficient  is  found  to  be, 

h  =  Ac  [  p  c  C  W  w -  Q  ] 
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From  the  definition  of  Nusselt  number  the  second  definition  of 
Nusselt  number  can  now  be  written  as 

(nu)2  -Bjia 

=  w_  _  F  w  (4.1) 

w9  w0*Re 

For  the  case  without  heat  sources  F  =  0,  equation  (4.1)  reduces  to 
equation  (3.8)  which  is 

(Nu)?  =  £ 

W0 

It  is  noted  that  the  accuracy  of  the  numerical  results  can  be  checked 
using  two  alternative  expressions  for  the  pressure-drop  and  heat 
transfer  parameters. 

The  effect  of  heat  sources  parameter  F  on  pressure-drop  and  heat 

transfer  parameters  is  demonstrated  in  Fig.  22  for  a  =  45°,  Pr  =  0.75, 

Re  =  100  and  Ra  =  40.  At  first  glance  the  pressure-drop  parameter  in 

the  form  of  fRe  appears  to  be  a  linearly  decreasing  function  of  the  heat- 

« 

generation  parameter  F  which  is  similar  to  that  reported  in  [7]  for  a 
different  definition  of  pressure-drop  parameter  for  combined  free  and 
forced  laminar  convection  in  vertical  rectangular  channels  with  internal 
heat  generation.  The  Nusselt  number  is  seen  to  decrease  appreciably 
with  the  increase  of  the  heat  sources  parameter  F.  The  difference  in  the 
relative  influence  of  heat  sources  parameter  on  the  parameters  fRe  and 
Nu  can  be  explained  from  the  velocity  and  temperature  profiles  shown  in  ‘ 
Fig.  19(a)  and  (b).  It  is  noted  that  at  F  =  41.5  the  heat  transfer  rate 
at  wall  is  zero,  and  with  the  further  increase  of  F,  the  Nusselt  number 
becomes  negative  indicating  reversed  heat  flow. 


, 
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A  comparison  of  the  heat  transfer  results  from  this  work  with  those 
obtained  by  a  perturbation  method  for  a  limiting  case  of  horizontal  tube 
[10]  is  shown  in  Figs.  23  and  24.  It  should  be  noted  that  the  Nusselt 
number  reported  in  [10]  is  based  on  the  mean  temperature  of  the  fluid 
with  (Nu)q  =  6  for  Ra  =  0  and  F  =  0;  whereas  the  Nusselt  number  from  this 
work  is  based  on  bulk  temperature.  It  is  expected  that  the  difference 
between  two  results  will  be  small  when  Ra  is  very  small.  However,  some 
difference  is  observed  for  F  =  10  in  Fig.  23,  and  Pr  =  5  in  Fig.  24.  In 
addition  to  the  difference  in  the  definition  of  Nusselt  number,  Fig.  23 
clearly  shows  the  typical  blow-up  trend  of  the  heat  transfer  results  with 
the  increase  of  Rayleigh  number  obtained  by  the  perturbation  method.  It 
can  be  expected  that  a  significant  difference  in  heat  transfer  predic¬ 
tions  exists  between  the  numerical  method  and  perturbation  method  when 
Ra  is  large. 

The  effect  of  tube  inclination  angle  on  flow  and  heat  transfer 
results  is  shown  in  Figs.  25  and  26  respectively,  with  heat  sources  func¬ 
tion  F  as  a  parameter  for  Pr  =  0.75,  Re  =  100  and  Ra  =  40.  A  close  exam¬ 
ination  of  the  effect  of  F  on  fRe  shown  in  Fig.  25  discloses  that  at  any 
angle  of  tube  inclination  the  value  of  fRe  is  not  a  linearly  decreasing 
function  of  F.  The  straight-line  relationship  shown  in  Fig.  22  is 
caused  by  a  rather  compressed  vertical  scale. 

The  effect  of  Rayleigh  number  on  pressure-drop  parameter  (fRe)  and 
Nusselt  number  for  a  =  0°,  45°,  and  90°  is  shown  in  Figs.  27  and  28 
respectively,  for  several  representative  values  of  F  with  Pr  =  0.75  and 
Re  =  100.  In  particular  applications,  the  high  Rayleigh  number  flow 
regime  is  of  interest,  and  the  effects  of  tube  inclination  and  heat 
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sources  are  both  seen  to  be  significant.  In  particular,  for  a  =  45°  and 
90°,  the  heat  sources  effect  is  predominant  when  the  Rayleigh  number  is 
large.  Fig.  28  shows  that  for  F  f  0,  the  effect  of  Rayleigh  number  on 
heat  transfer  in  horizontal  tubes  is  completely  different  from  that  in 
inclined  tubes  within  the  range  of  F  under  investigation. 

In  order  to  see  the  Prandtl  number  effect,  further  flow  and  heat 
transfer  results  for  a  =  0°,  45°,  and  90°  are  shown  in  Figs.  29  and  30 
respectively,  with  F  as  a  parameter  for  Pr  =  5  and  Re  =  100.  Note  that 
the  effect  of  Rayleigh  number  on  fRe  is  quite  small,  even  with  heat 
sources  for  horizontal  tube  a  =  0°with  Pr  =  5.  On  the  other  hand,  the 
effect  of  Rayleigh  number  on  Nu  is  significant  for  a  =  0°  at  Pr  =  5. 

This  observation  can  be  explained  from  the  fact  that  when  Prandtl  number 
is  large  and  the  secondary  flow  is  weak,  the  inertia  terms  in  the  momen¬ 
tum  equations  (2.7)  and  (2.9)  can  probably  be  neglected,  whereas  the 
convective  terms  in  the  energy  equation  (2.10)  are  important  because  of 
large  Prandtl  number  for  the  case  a  =  0°. 

It  should  be  noted  that  all  the  numerical  results  in  Chapter  III  and 
IV  are  obtained  by  using  a  mesh  size  of  M,N  =  28  and  a  relaxation  factor 
of  unity.  The  numerical  results  indicate  that  generally  the  first  four 
significant  figures  for  the  values  of  fRe  and  Nu  check  each  other  exactly 
for  the  two  alternative  expressions,  confirming  the  accuracy  of  the  num¬ 
erical  solution.  The  numerical  results  are  listed  in  Tables  1  to  11  for 
future  reference.  Tables  10  and  11  list  the  numerical  results  for  two 
alternative  expressions  for  fRe  and  Nu. 

The  computing  time  required  to  obtain  a  complete  solution  for  given 
values  of  the  parameters  Re,  Ra,  Pr,  F,  and  inclination  angle  a  ranges 
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from  2  to  3  minutes  up  to  Ra  =  400  on  the  IBM360/67  system.  The 
required  computing  time  increases  with  the  further  increase  of  Rayleigh 

o 

number,  and  generally  5  to  6  minutes  are  required  when  Ra  >  10  . 
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CHAPTER  V 

CONCLUDING  REHARKS 

1.  Based  on  the  present  numerical  results  it  is  now  evident  that  the 
perturbation  method  [1,  2,  10,  14]  as  used  in  the  literature  diver¬ 
ges  quickly  with  the  increase  of  the  characteristic  parameter.  The 
present  numerical  solution  indicates  that  for  any  combination  of  Ra, 
Re,  Pr,  and  F,  the  maximum  value  of  Nusselt  number  does  not  exist 
for  any  tube  inclination  angle  which  is  contrary  to  the  result  re¬ 
ported  in  [2,  14]  using  the  perturbation  method. 

2.  In  high  Rayleigh  number  regime,  the  tube  orientation  effect  on  flow 
and  heat  transfer  results  is  significant  in  the  neighborhood  of 

a  =  0°.  This  implies  that  at  high  Rayleigh  number  the  coupled  effects 
of  tube  orientation,  buoyancy  forces,  and  centrifugal  forces  [16] 
must  be  considered  for  laminar  forced  convection  heat  transfer  in 
helical  tubes. 

3.  The  numerical  solution  using  a  combination  of  boundary  vorticity 
method  and  line  iterative  relaxation  method  proves  to  be  very  effec¬ 
tive  for  combined  free  and  forced  laminar  convection  heat  transfer 

in  inclined  tubes  with  or  without  internal  heat  generation.  The  con¬ 
vergence  of  the  present  numerical  solution  is  ascertained  by  com¬ 
parison  with  the  known  numerical  results  for  the  two  limiting  cases 
of  horizontal  tube  (a  =  0°)  and  vertical  tube  (a  =  90°)  for  the  case 
without  heat  generation.  A  comparison  between  this  work  and  refer¬ 
ence  [5]  for  the  vertical  tube  with  uniform  internal  heat  sources  is 
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not  possible  because  of  the  different  definitions  for  heat-source 
parameter.  For  another  limiting  case  of  horizontal  tubes  with 
heat  sources  only  qualitative  comparison  with  reference  [14]  is 
possible  in  low  Rayleigh  number  regime  since  the  Nusselt  number  in 
[14]  is  based  on  mean  temperature  of  fluid  instead  of  bulk  temper¬ 
ature  used  in  the  present  investigation. 

4.  The  present  results  for  Pr  =  5  without  heat  sources  may  be  applied 
to  the  flat-plate  solar  collector  tubes  [15]  for  water  heating.  In 
view  of  the  rather  significant  tube  inclination  angle  effect  with 
the  increase  of  Rayleigh  number,  it  appears  that  the  determination 
of  the  suitable  film  heat- transfer  coefficients  in  solar  collectors 
based  on  theoretical  and  experimental  data  for  horizontal  tubes 
alone  is  inappropriate  and  unrealistic.  However,  the  determination 
of  the  film  heat-transfer  coefficients  in  inclined  solar  collector 
tubes  requires  the  solution  of  Graetz  problem  in  inclined  tubes. 

5.  The  effects  of  heat  sources  on  axial  velocity  and  temperature  pro¬ 
files  are  seen  to  be  considerable  in  inclined  tubes.  Reverse  flow 
can  occur  at  sufficiently  high  Rayleigh  number,  and  reversed  heat 
flow  appears  at  a  sufficiently  high  heat-source  parameter.  For 
given  values  of  the  parameters,  both  the  Nusselt  number  and  the 
pressure  drop  parameter  (fRe)  decrease  with  the  increasing  values 
of  heat-source  parameter  F.  For  a  given  value  of  F  (F  f  0),  the 
parameter  (fRe)  always  increases  with  increasing  values  of  Ra, 
whereas  the  Nusselt  number  generally  decreases  with  increasing  Ra 
in  high  Rayleigh  number  flow  regime  due  to  the  counteracting  effect 
between  heat  sources  and  impressed  heat  flux  at  the  tube  wall.  In 
the  case  without  heat  generation  (F  =  0),  both  flow  and  heat 
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transfer  parameters,  fRe  and  Nu,  increase  with  the  increasing 
value  of  Rayleigh  number. 

6.  The  effect  of  Prandtl  number  is  significant  in  low  Rayleigh  number 
regime.  The  Prandtl  number  effect  is  also  seen  to  be  significant 
when  the  inclination  angle  a  is  small.  For  any  given  value  of  the 
Prandtl  number,  an  asymptotic  behavior  for  flow  and  heat  transfer 
results  appears  in  high  Rayleigh  number  regime. 

7.  The  Reynolds  number  effect  appears  to  be  significant  in  the  neigh¬ 
borhood  of  horizontal  direction.  In  contrast,  for  the  vertical 
tubes,  the  flow  and  heat  transfer  results  are  independent  of 
Reynolds  number.  In  high  Rayleigh  number  regime,  an  asymptotic 
behavior  appears  for  flow  and  heat  transfer  results  for  given  values 
of  Prandtl  number,  heat  sources  parameter,  and  tube  inclination 
angle. 

8.  It  is  possible  to  apply  the  present  numerical  method  to  combined 
free  and  forced  laminar  convection  problems  where  the  body  forces 
are  caused  by  Coriolis  forces  or  buoyancy  forces  in  rotating  field 
with  a  body  force  orientation  effect.  Furthermore,  it  is  noted  that 
the  extension  of  the  numerical  technique  to  fully  developed  combined 
free  and  forced  laminar  convection  in  various  noncircular  inclined 
ducts  or  channels  is  obvious. 
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directions  <j>  =  0  and  tt. 
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Fig. 2(b)  Effect  of  Re  on  axial  velocity  profile  along  the 
directions  <j>  =  0  and  it. 
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Fig. 2(c)  Effect  of  inclination  angle  on  axial  velocity  profile 
along  the  directions  0=0  and  tt. 
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Fiq.3  Streamlines  and  isothermals  for  Pr  =  10,  Re  =  100, 
Ra  =  40,  and  a  =  45°. 
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Fig .4 (a)  Effect  of  Ra  on  Temperature  profile  along  the 
directions  <J>  =  0  and  tt. 
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Fig. 4(b)  Effect  of  Re  on  temperature  profile  along  the 
directions  d>  =  0  and  tt. 
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Fig.4(c)  Effect  of  inclination  angle  on  temperature  profile  along 
the  directions  <|)  =  0  and  it. 
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Fig. 5  fRe/(fRe)n  versus  Ra  with  a  as  a  parameter  for  Pr  =  0.75 

u  and  ReRa  =  4,000 
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Fig.  7'  Nu/(Nu)n  versus  Ra  with  a  as  a  parameter  for  Pr  =  0.75 

u  and  Re  Ra  =  4,000. 


Pr  =  0.75 
ReRa  =4,000 
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Fig. 8  Nu/(Nu)q .versus  Ra  sin  a  with  a  as  a  parameter  for  Pr 

0.75  and  ReRa  =  4,000. 
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against  perturbation  solution. 
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against  perturbation  solution. 
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Fis.13  fRe/(fRe)n  versus  Ra  with  Pr  as  a  parameter  for  ReRa 

=  4,000  and  a  =  45° 
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Fig. 15  fRe/(fRe)n  versus  a  with  Ra  as  a  parameter  for  Pr  =  0.75 

and  ReRa  =  10,000 


Fig. 16  Nu/(Nu)n  versus  a  with  Ra  as  a  parameter  for  Pr  =  0.75 

u  and  ReRa  =  10,000. 
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Fig. 17  Comparison  of  heat  transfer  results  from  this  work  with 
those  from  perturbation  method  [2]. 
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from  perturbation  analysis. 
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Fig.  19(a)  Typical  velocity  profiles  with  F  as  a  parameter 
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Fig.  19(b)  Typical  temperature  profiles  with  F  as  a  parameter 


64 


0=0  -= -  j  - *-  11=  0 


S- 

<d 

+-» 

<u 

s 

nj 

S- 

nj 

Q. 

nj 


Fig.  20(a)  Typical  velocity  profiles  with  Ra  as 


Pr=  0.75 
Re=  100 
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Fig.  20(b)  Typical  temperature  profile  with  Ra  as  a  parameter 
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Fig.2l(a)  b]stnbution  of  secondary  velocity  component  u  along  <p  =  0  and 
with  Ra  as  a  parameter  y 
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Fig .2l(b)  Distribution  of  secondary  velocity  component  v  alonq  d>  =  w? 

with  Ra  as  a  parameter  ■  Y  ' 
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Fig. 22  fRe/(fRe)0  or  Hu/(Nu)Q  versus  F. 


PERTURBATION  ANALYSIS 
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Fig. 24  Comparison  of  heat  transfer  results  between  numerical  solution 
and  perturbation  solution. 
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Fig. 25  fRe/(fRe)n  versus  a  with  F  as  a  parameter 
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Fig. 26  Nusselt  number  ratio  versus  a  with  F  as  a  parameter 
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Fig. 2 7  fRe/(fRe)Q  versus  Ra  with  F  as  a  parameter  for  Pr  =  0.75 

and  a  =  0°,  450  and  90° 
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Fig. 2 8  Nusselt  number  ratio  versus  Rayleigh  number  with  F  as  a  parameter 
for  Pr  =  0.75  and  a  =  0°,  45°  and  90° 
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o 
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Fig.  29  fRe/(fRe)n  versus  a  with  F  as  a  parameter  for  Pr  =  5  and 

u  a  =  0°,  45°  and  90° 
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Fig. 30  Nusselt  number  ratio  versus  Rayleigh  number  with  F  as 
parameter  for  Pr  =  5  and  a  =  0°,  45°  and  90° 
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APPENDIX  A 

REDUCTION  OF  GOVERNING  DIFFERENTIAL 

EQUATIONS  TO  DIMENSIONLESS  FORM 


The  governing  conservation  equations  for  mass  momentum  and  energy 
in  cylindrical  coordinates  for  a  steady  fully  developed  upward  laminar 
flow  in  an  inclined  tube  (see  Fig.  1)  using  the  simplifying  assumptions 
made  in  Chapter  II  are: 

Continuity  equations: 


|r  (RU)  +  =  0  (AJ) 

Momentum  equations: 

R-direction, 


P  (U 


V  3U 
R  3<J> 


-  p  g  cos  a  cos  (J)  - 


3P 

3R 


,  /32U  ,  1 

+  ]i  ( — w  +  - 


3R 


R 


3U 

3R 


1  32U 

R2  3(J)2 


U__  2_  3V\ 

R2  "  R2 


(A. 2) 


(^-direction, 


P  (u 


8V  ,  V  3V  .  UV V 
3R  R  3<p  R  1 


1  3P 

=  p  g  cos  a  sin  <t>  -  ^ 


2 

/3  V  , 
+  p  (— o  + 

3R 


1  3V  1 
R  3R  r2 


2_  3U 
R2  8* 


(A. 3) 


Z-di recti  on, 


P  (U 


3W 

3R 


V  3W\ 
R  3<jr 


=  -  p  g  sin  a  - 


fr*^ 


32W  ,  1  9W  . 

3R2  R9R  R2 


2 

3  W' 

2' 

3<j) 


(A. 4) 


1 
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Energy  equation: 


r  #„3T  ,  V  9T 
p  Cp  U3R  R  3<J) 


W§I)  = 


(£+1 

3r2  R  9R 


3T  ,  1 


A  +  A)  +  Q  + 


$ 


R  3<J> 


3Z 


where  the  viscous  dissipation  function  $  is  defined  as 


(A. 5) 


A  o  / / 9U \ 2  rl 
$  =  2y  j(^)  +  [- 


R 


<$+w2 


u 


fi  im.\2 

lR  dp' 


+  (3W)2  +  rl  3U  +  R  d_  /Vn2l 
3R  LR  34*  +  R  3R  (R,J 


(A. 6) 


and  Q  is  the  rate  of  internal  heat  generation  per  unit  volume  which  will 
be  considered  to  be  constant  in  this  study.  In  order  to  normalize  the 
governing  equations,  it  is  convenient  to  introduce  the  following  charac¬ 
teristic  quantities: 

U  =  Ucu,  V  =  Vcv,  W  =  Wcw,  R  =  ar,  Tw  -  T  =  0C0  (A. 7) 

From  the  continuity  equation,  the  secondary  velocity  components  U  and  V 
can  be  considered  as  the  same  order  of  magnitude  as  indicated  below. 


Setting,  V  =  U  ,  the  continuity  equation  is  reduced  to 
c  c 


a(ru)  ,  ay. 

3r  3<f> 

Using  Boussinesq  approximation 
term  can  be  expressed  as, 


=  0  (A. 8) 

the  density  variation  in  buoyancy 


P  =  Pw  [H3(TW-T>]  =  pw  [1+3  0c0]  (A*9) 

It  is  recalled  that  the  axial  pressure  gradient  is  assumed  to  be  constant 
and  the  axial  temperature  gradient  is  constant  due  to  the  imposed  thermal 


. 
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boundary  conditions,  namely 


gp 

■  (gY  +  Pw  9  sin  a)  =  A  =  constant 


£L 

az 


=  c  =  constant 


The  characteristic  axial  velocity  W  can  be  taken  as  the  maximum  velo- 

V 

city  in  the  pure  forced  convection  case,  that  is, 

W  =  W  =  ^  =  2  W 
c  max  4 y 

and  Reynolds  number  Re  will  be  defined  using  mean  velocity  as, 


Re  = 


vJ  2a  A  a3  Wc  a 


V 


4pv‘ 


By  introducing  the  above  characteristic  quantities,  the  axial  momentum 
equation  (A. 4)  can  now  be  written  as 


V2w  =  (u|^  +  -  |t)  +  [9  6  M'C  3  ]  6  sin  a  -  " 

L  v  J  V“gr  r  L  V  W_  J  C 


y  W 


(A. 10) 


c  c 

It  is  possible  to  set  the  coefficient  of  the  convective  terms  to  be 

order  one,  and  we  have 

u 

c  a 

It  is  noted  that  each  term  of  the  equation  (A. 10)  is  seen  to  be  of  the 

same  order  of  magnitude. 

Similarly,  equation  (A. 5)  can  be  written  as 

2 


+  <P 


(A. 11) 


a  ■  <$**!!>-  ^ 

The  conduction  terms  and  the  convection  term  in  the  energy  equation  may 

be  considered  to  be  the  same  order  of  magnitude  leading  to. 

p  C  U  a  y  C 
H  p  c  _  _ J)  =  pr 
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W  C  p  C  a' 
-JL .. .  P  . 

k  0. 


-  1  >  0„  =  Re  C  a  Pr 

c 


The  dimensionless  viscous  dissipation  function  is 

2 


,  _  a  <£> 


= 


k  0 


(A. 12) 


_  r  f  2  r  / 5 u \  2  1  / 3 v  ,  \2-i  ,  n  r  / 1  3W\2 

-  Ec  \Re  [(3?)  +  T  (3?  +  u)  ]  +  Re  a*) 


.  /3wn 2-i  J r  1  3u  ,  3  /V\-i2l 

+  W  ]  +  r¥  c7  3*  +  *3?  (?)]  j 


where  = 


U 


c  cp  C  a 


=  Eckert  number. 


If  E  «  1,  the  viscous  dissipation  terms  in  equation  (A. 11)  can  be  neg- 
c 

lected.  The  dimensionless  heat-source  parameter  is  defined  as, 


F  = 


Q  a 


pCpvC 


Then  the  normalized  energy  equation  becomes, 


v2e  =  Pr  (u|£ 

Substituting  equation  (A. 12) 
equation  becomes 


+  *!!)  -w  +  F 

r  3(f)7 


Re 


into  equation  (A. 10) , 


(A. 13) 

the  axial  momentum 


2 

V  w  = 


u 


3w 

3r 


+  * 


3w 

3<j) 


+  Ra  0  sin  a  -  4 


(A. 14) 


where  Ra  =  3  9  C  a^/vK,  is  Rayleigh  number. 

The  pressure  terms  in  the  equations  (A. 2)  and  (A. 3)  can  be  elimin¬ 
ated  by  differentiating  equation  (A. 2)  with  respect  to  <p  and  equation  (A. 3) 
with  respect  to  R,  and  then  reducing  them  to  one  equation.  Introducing 
the  dimensionless  stream  function  and  uhe  vorticity  as, 
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u  =  -M- 

r3<J)  ’ 


=  M. 

3r 


5  =  V> 
2 


(A. 15) 


where  V*"  =  — +  1  — +  V  - 

where  V  ar  +  r  9r  +  ?  . .2 

r  dtp 


the  momentum  equations  (A. 2)  and  (A. 3)  are  reduced  to: 


V2£  -  u||-  +  ~  |~  +  ReRa  (|~  sin  <J>  +  ~  ~  cos  <j>)  cos  a  (A.  16) 

A  set  of  partial  differential  equations  (A. 13)  to  (A. 16)  is  the  normal¬ 
ized  equations  which  have  been  solved  numerically  in  this  study.  For  the 
problem  under  consideration,  five  independent  parameters  Re,  Ra,  Pr,  a 
and  F  appear  in  the  governing  dimensionless  differential  equations. 
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APPENDIX  B 
FORTRAN  PROGRAM 


»*'  IN  CCL •  6:  THIS  LINE  IS  A  PART  CF  PREVIOUS  STATEMENT 


C 

c 

c 

c 


THIS  IS  A  PROGRAM  FOR  SOLVING  COMBINED  FREE  AND  FORCED 
CONVECTION  WITH  HEAT  GENERATION  IN  INCLINED  TUDES 


c 

EXPLANATION  OF  SYMBOLS 

IN  THIS 

PROGRAM 

c 

c 

FN 

• 

FEAT  GENERATION 

PAR AMENT  NUMBER 

c 

M 

• 

* 

NO.  OF  DIVISION 

IN  THE 

R-D I RECTION 

c 

N 

• 

• 

NO.  CF  DIVISION 

IN  THE 

CITA  DIRECT 

c 

PR 

• 

• 

PRANDTL  NUMBER 

c 

RA 

♦ 

• 

RAYLEIGH  NUMBER 

c 

RE 

• 

-• 

REYNOLDS  NUMBER 

c 

S 

• 

* 

STREAM  FUNCTION 

c 

VC 

• 

♦ 

VELOCITY  IN  THE 

CITA  DIRECTION 

c 

VQ 

* 

• 

VGTIC I TY  FUNCTION 

c 

VR 

• 

• 

VELOCITY  IN  THE 

RAD  I AL 

DIRECTION 

c 

VX 

* 

-* 

VELOCITY  IN  THE 

AXIAL 

D  IRECT ION 

C 

c 


C 

c 


DIMENSION  VX ( 3 1  *3  1)  »  T  {  3  1  *31)  ,  VR ( 3 1  *  3 1  )  »VC(J1  *31  )  ,BVTS( 

£  31*31). 

1AVX(21*31  )  *  C  V  X ( 3  1  *31)  .  P  V  X ( 3 1  .31  )  »S(31  *31  )  *VL  (31  .31  )  ,AT 

*  (31*31). 

3Cy(3i  ,31)  *PT(31  *31)  *  A  S ( 3 1  *31  )  »CS{31  *31  )  *  P  S ( 3 1  *  31  )*TT(3 

CCITA  (31  )»R (31  )*ZA(3  1  )»ZB(31 ) *ZC(31 )  »VVX( 31 *31 )  *  SS ( 3 1 ,3 

❖  1  ) 

INTEGERS  M.NiMI  *  N  I  *  N  D  *  M  D  *  M  S  TD  *  NO 

COMMON  HR ,HR2.HC,HC2. AF , R AC , R A , PR  * S TD *  OK , ZA , ZB . ZC * R . C I 
*  T  A  ,  E  V  T  S  . 

1VX,T,VR*VC.CMEC-1,GMES 

CUM  MON  M  *  N  *  Ml  I  *  N  1  *  N'D  *  MD  *  MSTD  *  NO 
99  READ ( 5  »  ICC  *  E  N  D  —  1 50)  M.N.STD.CK 

REA  C  IN  THE  GIVEN  VALUE 
READ!  5  *  2  C  1  )  RE » R A  ,  AF . PR  * FN 


RAC-RE*RA 


' 
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c 

c 

c 


F  .V  =  F  N  /  R  E 

RE£D  (  3  ,  250  )  C.V.EG1  ,  C  NE  S  •  N*S  T  O  ,  MO 
P 1=3  •  14  1  3  927 
AG  =  AF 


AF  =  P  I  *AF/ 1  EC 


33  MD=W 

- 1 

V  I  =  V  +  l 

ND  =  N 

- 1 

N  I  .=  N  +  1 

N0  =  0 

NCC  = 

5  C 

MO=^-C  +  1 

I  F  {  M 

,GT  . 

2  1 

•  OR .^C  ,  G  T 

DO  10  J= 

1 

f 

NI 

DO  10  1= 

1 

9 

N  I 

VX  (  I 

1  J  )- 

0 

• 

0 

V  VX  ( 

I  ,  J  ) 

- 

0 

.0 

T  (  I 

»  J  )  = 

0 

• 

C 

TT  (  I 

,  J)  = 

0 

• 

0 

Si  l 

»  J  )  - 

c 

• 

c 

SS  (  I 

♦  J  )  = 

0 

• 

0 

VO  (  I 

,  J)  = 

0 

• 

0 

VR  (  I 

,  J  )  = 

0 

• 

0 

VC  (  I 

»  J  )  — 

0 

-a 

0 

10  CGMIME 


1  )  GO  TO  3  4 


READ  IN  TEE  INITIAL  VALUES 


34 


9 


READ (3,200) 
READ ( 3,200) 
READ { 5,200) 
READ  (3,200 
READ ( 3 ,200 ) 
RE  ADI  3 , 2  C  C  ) 
READ ( 3 , 2  C  C ) 
READ ( 5 , 2C0 ) 
READ! 3 . 2  0  C ) 
ER= 1  .  C/N 
HR2=ER**2 
F=F.v*HR2 
E  C  =  P  I  /  N 
HC2=hC**2 


(  ( VX<  I  ♦ J)  , 
(  {  V  V  X  (  I  ,  J  ) 
{  {  T  (  I  ,  J  )  , 
(  (TT {  I , J  )  , 
(  <  S  {  I  ,  J  )  , 
((SS(I.J), 
(  ( VO (  I  , J )  , 
UVR(I.J), 
(  ( VC (  I  ,  J )  , 


1  =  1  ,MI  )  , J  = 
•  I  =  1  •  M  I  )  ,  J 
[  =  1  , M I  )  »  J  = 
I  =  1  ,  M  I  )  ,  J  = 
1=1 ,MI ) , J= 
1  =  1  ,MI  )  » J  = 
I  —  1  ,  M  I  )  ,  J  = 
I  =  1  , M I  )  ,  J  = 
I  =  1  ♦ M I  )  , J  = 


C  I  T  A  (  1  )  =  0 . 0 
DO  9  J  =  2  ,  N  I 
C  l  TA ( J  )  =  C I TA ( J-l  ) +  EC 


R  (  l  )  =  0 .0 


DO  3  1=2,  M 


R(  I  }  =  R  (  I— l  )  +HR 
ZA (  I )  =  E R / ( 2*R (  I  )  ) 

ZG(  I  )-(ER/(EC*R(  I  ))  )  *  *  2 
ZCC I )-ER**2/(2*R(I)*HC) 


8  CONTINUE 


»NI  ) 

1  ,NI  ) 

•  N  I  ) 
,N  I  ) 

,  N  I  ) 

,  N  I  ) 

.  N  I  ) 

,  N  I  ) 

•  NI  ) 
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DO  7  J  = 1 ,NI 
DO  7  I  =  2,. VI 

evrs( i »j)=( i+ze(  i )) *(— 2) 

AS (  I  , J )- 1 -ZA {  I  ) 

CS {  I i J )= 1 +ZA {  I  ) 

7  CONTINUE 

IF  {  f/  .L  T  ,2  1  )  GC  TC  1 
CO  2  2  J-1,NII,2 
DO  2  3  I=2*jV,2 

VX{  I  , J)=0.5*(VX{  1-1  *  J  )  +  VX  {  I  4-  1  ,  J  )  ) 

VVX(  [  ,  J )-0*5*(VVX(  1-1  »  J  )  +  V  V  X  (  14-1  ,J)  ) 

T(I,J)=0.5*{  T (I— I»J)+  T(I+1,J)) 

TT< I ,J)-0*5*(TTiI-l,J)+TT(l4-l,J)) 

S(I,J)=0.5*(  S {  I - 1 , J ) +  S (  I  4-  1  , J ) ) 

SS(  I*J)=0.5*(SSCI-1,J)+SS(  I  +  1,J)) 

VO(  I  .J)  —  0»5*(VG(  I  —  1  »J)  +VC{  1  +  1  »J)) 

V  R  (  l*J)=0.5*{VR(  I - i *  J ) +  V  R (  I  +  1  *  J ) ) 

VC(I  *J)=C.5*(VC<I-1  ♦ J  )  +  V  C (  14-1  ,  J  )  ) 

23  CONTINUE 

DO  24  1=1  ,  JV  I 

DO  24  J  =  2  » N  »  2 

VX{  I  ,J)=0.5*(VX{  I  ,  J-l  )  4-VX(  I  ,  J4-  l  )  ) 

V  V  X  {  I  ,J)=0.5*(VVX(I,J-1)4VVX{I,J4-1)) 

T(I,J)=0.5*(  T(I,J-n+  T  {  I  *  J+  1  )  ) 

TT<  I  »J)  =  0«5*{TT {  l,J-l  )+TT(  I  , J+ 1  )  ) 

S(I,J)=0.5*(  SUiJ-1  )f  S(I,J+1)) 

SS{  I#J)  =  C*S*(SS(  I  »  J-  l  ) +SS ( I  ,  J4-  1  ) ) 

VO(  I»J)=0.5*{VG(  I*J-l)+VO(I*J+l)) 

VR  <  I  *J)-0.S*(VR{  I  ,  J-l  )  +  V  R  (  I  ♦  J+  1  )  ) 

VC{  I , J  )  =  0 .5* ( VC (  I  . J- 1  ) 4 VC (  I  , J+ 1 ) ) 

24  CONTINUE 
1  I  M=  1 

C 

C  CO  LOOP  OF  25  IS  TO  CALCULATE  THE  COEFFICIENT  OF  LINE 

C  ITERATION 

C 

DO  25  J=l,Nl 
A  S ( M  I  »  J  )  =  2  •  0 

AVX { 2  ,  J  )= 1-ZA ( 2  )  +  0.5*HR*VR (2* J) 

CVX(2,J)-l4-ZA{2)-HR*VR(2,J)/2 

PVX(2,J)=CVX(2»J)/0VTS(2»J) 

AT{2,J)=1-ZA(2)4C*5*PR*HR*VR(2,J) 

CT<2,J)-1+ZA<2)-RR*PR*VR(2.J)*0*5 

PT(  2  »  J  )  =  C  T (2  *J)/QVTS(2»J) 

PS(2,J)=CS<2,J)/EVTS(2,J) 

DO  25  1-3,  JV 

AVX {  I  »J  )  =  1  — Z  A  (  I  }4-0.5*FR*VR<I  ,J) 

CVX (  I  ,  J  )=  1  4-ZA  (  I  )-0  •  5*Ff«*VR( I  *J) 

P  VX  (  I  ,  J  )  •=  C  V  X  (  l  ♦  J  )  /  (  0  V  T  S  (  I  ,  J  )  -  A  VX  {  I  ,  J  )  *  P  V  X  (  I  -  1  ,  J  )  ) 

AT (  I  ,  J  )=  1-ZA (  I  )+C*5#HR$PR*VR (I  »  J ) 

C  T (  I  ,  J  )  =  l  4  Z  A  (  I  )-0.5*FR*PR*VR(I»J) 

P  T {  I  , J )-CT (IfJ)Z(BVTS(I  ,  J ) -AT (  I  ,  J ) * P T (  I - 1  »  J )  ) 


* 

■ 


PS (  I  .  J  )  =  CS (  I.J)/(8VTS(I,J)-AS(  l  ,  J  )  *PS (  l  -  1  , J )  ) 

25  CONTINUE 

CALL  SUE V> ( V VX  ♦  A V X  ,  C V X » F V X ) 

CALL  SUET < TT ,AT  ,CT  ,PT.F ) 

CALL  SUEVS(S.VC.$S,AS,CS,PS»AVX,CVX,PVX) 

CALL  SUEVRC ( S  .ERROR ) 

IF  {  ERRCR  .LT  .  CK  .  CR  ,N’C  •  GT  .  VSTO  )  GG  TO  60 
N  0  — NC  f  1 

I F  (  N  C  « N E  * N C G  )  GG  TO  1 
WRITE (6*101)  EFFCH 
WRITECC.llC)  N  G 
NCC-FCC  +  50 
GO  TC  1 

60  CALL  TANK(SVXP»SVXT, SGT, SGVX. BUKT.FREi  .FRE2, UNI  .UN2.RE 
*  RA  .  RFE  1  ,  V- 

CRFE  2  ,  RNU  1  .RNU  2  .  ARNU *S  TF ,FRE J , A VEN , AVEF ,  RFE l 2  *  RFE3 , FM ) 


IF( M .LT 

20)  GG  TG  5 

W  R  I  T  E  1  6 

105) 

WRITE (6 

113) 

<  ( VX{  I , J )  *  1=  l  *  M  )  . J  = 

WP  ITE (6 

ICC) 

WRITE  (6 

113) 

(  (T (  I  *J)  •  I=1*M)»J=1 

W  R  I  T  E  ( e 

1  C  7  ) 

U  R  I  T  E  (  € 

113) 

(  (  V  0  (  I  *  J  )  .  I  =  2  »  M  1  )  .  J 

WRITE (6 

1  C€  ) 

W  R  I  T  E  (  e 

113) 

(  ( S (  I  , J  )  ,  1  =  2  ,MI  5  »  J  = 

W  R I  T  6  (  6 

109) 

w  R  IT  E  (  6 

113) 

(  (VR( I  * J) *  1-1  .  M)  .  J  = 

W  R  I  T  E  {  6 

l  l  0  ) 

W  R  I  T  E  (  6 

113) 

(  ( VC  <  I  »  J )  .  1= 1 ♦ M  )  »  J  = 

IM=5 

W  R  I  T  6  {  t 

202  ) 

N  ,  N 

W  R  I  T  E  <  C 

203  ) 

STD  *GK 

W  R  l  T  E  (  e 

2  04  ) 

RAC  .RA.PR.RE.FN 

WRITE (6 

104) 

AG 

WR I TE  (  6 

10  1) 

ERROR 

W  R  I  T  E  (  6 

116) 

NC 

w  r  i  r  e  ( e 

2  C  5  ) 

CMEC1  ♦ GME  S 

WR I TE  (  6 

2  10) 

WRITE  (6 

111) 

□  UK  T 

WRITE  (6 

212  ) 

STP 

WR  I  TE  (  e 

112) 

SVXF 

W  R IT  E  (  6 

114) 

ERE  1  »FRE2 

WR  l  TE  (  6 

12  1) 

RFE1  » RFE  2 

W  R  I  T  E  (  6 

2  15) 

RFE  12 

WR  IT  E  (  6 

2  14) 

FRE3 . RFE3 

WR I TE ( e 

115) 

UNI  .  U  N  2 

W  R  I  T  E  (  t 

2  12) 

A VEN. AVEF 

W  R  I  T  E  (  C 

1  22  ) 

RNU  1  . RNU2 »  ARNU 

W  R  I  T  E  (  C 

120) 

SVXT  ,  SGT , SGVX 

Vi  R  I  T  E  (  6 

117) 

RERA 

WRITE  (6 

2  11) 

IF  (  MC *LT •5*0R«M#GT«21  )  GC  TO  75 

. 


, 
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c 

c 


c 

c 

c 


V?R  I  TE  (  7 .200  ) 
toR I TE ( 7  *  20  C  ) 
*RITE(7,20C) 
vvRITE  (7  ♦  2  0  C  ) 
*RITE(7»20C> 
W  R  I T  E  (  7 . 2  0  0  ) 
WR  I  TE  (  7 .200  ) 
V.RI  TEC  7  »2CC) 
W  R  1  T  E  (  7 , 2  0  0  > 


(  ( V  > (  I  i J )  «  I = 1  .MI)  *  J  =  1 
(  (VVX(  l  »  J)  .  1=  1  )  »J  = 
( (  T ( I . J ) ♦ 1= l ,M I ) . J=1 
((TT(  I i J )  ♦  I- 1 »MI  ) , J= 1 
( (  S  (  I  .  J  )  •  1=  1  .M  I  )  .  J=1 
( ( SS(  I  * J )  i 1-1  *MI )  » J=  1 
(  (VC(  I»J)»  I  = 1 »  W  I  ) , J= 1 
(  (  V  R  (  I  *J)  .1  =  1  .MI)  »  J=  1 
(  (VC(  I»J)  *  1  =  1  .MI)  »J=1 


I  F  (  M  •  G  T #21  )  GC  TC  4  0 


DO  El  J=  1  *NI 

L  =  2*  (  N  I- J  )X+  1 

J  J  =  N l  +  1- J 

DO  3  1  1  =  1  .MI 

K -2*  ( MI  -  l  >  + l 

I  l  =  M  I  +  1  -  I 

VX( K  »  L  )  =  V  X (  I  I  ,  JJ  ) 

V  V  X  (  K  ♦  L  )-VVX (  I  I  ,  J J  ) 
T ( K  »  L ) -  T(II.JJ) 

TT (K  ,L  )  =  TT (  I  I  .  J J  ) 

S (K  ,L  )  =  S  (  I  I  ,  J  J  > 

SS ( K .L )  =  SS (  I  I  ♦  J  J  ) 

VO ( K , L  )  — VC (  I  I  . J J  ) 
VR{K,L)  = VP <  I  I  .  J J  ) 

VC ( K . L )  =  V C  (  I  I*JJ) 

2  1  CONTINUE 


ENLARGE  TEE  MESE  SIZE 


M  =  2*N 
N  =  2*N 
GO  TG  22 

40  DU  4  1  J= 1 . N  I  •  2 
L  = {  J+  1  )/2 

DO  41  1=1. MI. 2 

K  =  (  I  +  1  )  / 2 

VX (K ,L ) - VX (  I  .  J  ) 

V VX ( K  ,L  )  =  V VX (  l.J) 

T ( K . L  )  =  T(I.J) 

TT(K,L)~TT(I .J) 

S ( K  .  L  )  =  S(I.J) 

SSI K  .L  )  =  SS  I  I  . J ) 

VQ( K  .L  >  = VC  < I  *  J ) 

VR (K.L  )=VR (  I  *  J  ) 

VC(K»L)=VC( I »J) 

41  CONTINUE 
GO  TC  <SS 

100  FORM  AT ( 2  I  4  ,2F  l  0 *6 > 

10  1  FORMAT!  'O'  ,50X,< ERROR- 15X,bl4.6> 


.  N  I  ) 

1  ,N  I  ) 

♦  N  I  ) 

•  N  I  ) 

»  N  I  ) 

» N  I  ) 
iNI  ) 

»  N  l  ) 

»  N  I  ) 
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104  FOR?//*  T  (  ’  C  '  .  4  OX  .  *  A  NGLE  IN  DEGREE' 


105  FORMAT (  *  1  *  «  • 

*  1  1  C  N  '  ) 

107  FORMAT!  *1*  ,  ' 

1 06  FORMAT !  »  1  »  ,  » 
1  08  FORMAT (  »  1  '  ,  • 
1  C9  FORMAT (  '  1  '  ,  • 
110  FORMAT!  •  1*  ,• 

*  ) 


AXIAL  VELOCITY'  *  3  X  . 

VOT  IC  I  TY  »  ) 

TEMPERATURE  *  ) 

STREAM  FUNCTION') 
SECCNOARY  VELOCITY 
SECONDARY  VELOCITY 


,F 1 0 . 1  ) 

'  I  FROM  1  TO  M  J  FROM 


IN  R  —  DIRECT  ICN  VR '  ) 

IN  C  ITA-DI RECT ION  VC* 


111  FORMAT! * C ', 1CX, »0ULK  TEMPERATURE... 

^  ...  ,» «  .  * 

C .F 15.6) 


112  FORMA T  !* 0  ',  10X ,  •  AVERAGE  VELOCITY  IN  AXIAL  DIRECTION... 

*  . . 

C  ♦  F  1  5 . 6  ) 

1  13  FORMAT !  *  * , 7E  18. 7 ) 

114  FORMAT  {'O'  .lOX.'FREl  FRE 2.. . 

^  ... 

C.F15.6.F25.6) 

115  FORMAT  !  » C »  ,  l OX ,  • NUSSELT  NUMBER  NU  1  NU  2 • »  •  . . 

...  .  .  ' 

C  »F 15 .6  ,F25  .6  ) 

116  FORMAT (* 0 *. 45X ,' NUMBER  OF  ITERATIONS*  #  1  OX « I  l 0 ) 

117  FORMAT {  • C *  ,  1  OX  .  • PER  A .  .  .  ...  . . ,  ... 

^  ...  ...* 

C  .  F  1  5 . 6  ) 


120  FORMAT!  » C '.  1  OX SVXT  GRADIENT  OF  T  AND  V  AT  WALL. 

*  . 

C .FI  5.6 »  2  F  2  5  •  6  ) 

121  FORMAT (' 0 '.  1  OX  ,  »FRE  1/FREC  FRE2/FRE0 . 

4  ...  ...' 

C.F1S.6.F25.6) 

122  FORMAT {  ' C '  .  I CX  ,  ' NU1 /NUO  NU2/NU0  AVERAGE  OF  THEM 

*  . • 


C  *F1  5 .6 . 2  F  2  5  •  6  ) 

200  FORMAT! 20 A4) 

2  0  1  FO  RM  AT ! 3  F 1 0 . 1  .F1C.4.F1C.4) 

202  FORMAT!'  1*  »/////» 4 OX.  'MESH  S  I  ZE *  ♦  I  5 » 2 X  .  ' 0 Y •  ,  I  5 ) 

203  FORMAT!  'O'  .  3  0  X  ,  'ERROR  L I M  I  TE C '  , 2E2 0 . 6 ) 

204  FORMAT!  '0*  .25X.  'RAC'  .F1C.1.5X,  'RA'  .F10.  1  »SX,  'PR*  .F10.4 


*  , 5X  ,  *  RE  »  . 

CF10.l*SX.'F».F10.4) 

205  FORMAT! 'O'. 40X. 'REDUCING  FACTOR  IN  V  AND  T'»5X»F5.2.5X 

*  .'IN  S  *  .  SX  ♦ 

CF5 .2  ) 

210  FORMAT! *0* ,///* it**************************** y*  ******** 

*  6'  *  FEAT 

CTRANCFLR  RESULT  *************************£*'******* 

*  ****♦•*///) 

211  FORMAT!  *C'  ♦///*•**********  **************************** 

*  $  *  4  *  *  $  * 

CENp  a******************************************** 

*  *  »  , 


' 

■ 

' 


U)  U) 
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212  FORMAT (' C" t 10X •• MEAN  TEMPERATURE . 

^  •••* 

C.F15.6) 

213  FORMAT  (* C  •*  10X ,» AVERAGE  RUSSEL T  NC  AND  FRICTION  FACTOR 

*  . . * 

C.F  15  .e ,F25.C) 

2  14  FORM AT  (  •  C *  ,  1  OX ,  » PRESSURE  DROP  AND  IT  RATIO. . 

C.F15.C.F25.E) 

2  15  FORMAT  (•  C  •»  10X*  •  AVERAGE  GF  FRE1/FREO  AND  FRE2/FRE0. •  •  • 

*  . •••* 

C  »F 1 5  .  e  ,F25 .6 ) 

0  FORMAT ( 2F5 .2 , 2  I  5 ) 

0  S  TCP 
END 
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C  SUBROUTINE  SUEVX  IS  TG  SOLVE  THE  MOMENTUM  EGUATIGN  OF 

C  (2-9)  IN  CHAPTER  2 

C 

c 

c 

SUBROUTINE  SUEVX  (MX  *AVX  ,  C  V  X  ,  P  V  X  ) 

DIMENSION  V  X  (  2  1  »  3  .1  )  »  T  {  3  1  »  3  1  )  *VR(31»31)  ,  VC  (  3  I  •  3  1  )  *  8  V  T  S  ( 

*  31,31), 

1  A  V  X  (  2  1  ,  2  1  )  »  C  V  X  (  2  1  ,31)  , F  V  X ( 3 1  ,31)  , G  V  X ( 3 1  ,31  )  , D  V  X ( 3 1  ,31  > 

*  , 

CCITA(31),R(31  )  ,  Z  A ( 3  1  )  ,ZE(3l),ZC(31),VVX(21,3l) 

I  NT EGER *4  M ,N ,  MI  ,  N  I  ,ND,MD,MSTD,NO 

COMMON  HR , HR 2 , HC ,FC2, AF,RAC , RA , PR  *  STO,CK , ZA  ,ZB ,ZC,R ,C I 

*  TA,BVTS, 

1  V X ,  T , VR  ,VC  *  0  M  £  G  1  ,CMES 
COMMON  M  ,  N , M I , NI  ,  ND, MD , MS TO , NO 

VVX (1,1  )  =  0 .25* {  ( 1  —  H  R  *  V  R  (1,1  3/2) *VX (2,1  )  +(  1 +HR*VR(  1,1)/ 

*  2 ) *  V  X ( 2 , NI  ) 

CH2*VX(2,N/2+l)+(4-RA*T(  1  ,1  )  *SIN(  AF  )  )*HR2) 

DO  9  J  =  1  ,N  I 

VVX (  1  , J  )  =V VX{  1,1) 

IF ( J,EC  .  1  .OR  .J  .EG  .N  I  )  GO  TO  6 

DVX( 2,J)  =  (RA*T ( 2 , J ) *  S I N { A  F )  —  4 ) *  HR  2  — ( Z  B ( 2 )  +  Z  C ( 2 ) *  VC ( 2 , J 

*  )  ) * VVX ( 2  ,  J-  1 

C)-(ZE(2)-ZC(2)*VC(2,J) )*VVX(2,J+l)-AVX(2,J)*VVX( 1 , J ) 

6  IF(J.EQ.l)  DVX ( 2  ,  J )  =  ( RA*T ( 2, J ) *SI N ( AF ) -4 ) *HR2-2*ZB (2 ) * 

*  V VX  ( 2  ,  J  +  1  )- 
CVVX(  1  ,  J)*AVX(2,J) 

IF(J.EG.NI  )  DVX(2,J)-(RA*T(2,J)*SIN(AF)— 4)*HR2  —  2*ZB(2) 

*  *VVX ( 2 , J  —  1  ) 

C-VVX(  1  ,J)*AVX(2,J) 

GVX ( 2  ,  J  )-CVX ( 2,J)/BVTS(2,J) 

DO  10  1=3, M 

IF(  J.EC.l  .GR.J.EC.Nl  )  GO  TO  7 

CVX(  1  , J )  =  ( RA*S  IN ( AF ) *T(  I  , J )-4 ) *HR 2- ( ZB  <  I ) +ZC (  I ) *  VC (  I  , J 

*  )  ) *  V VX  (  I  ,  J- l 

C )  —  ( Z  E  (  I  )  —  ZC (  I  )  *  V  C (  I  ,J))*VVX(  I  »  j  +1 ) 

7  IF (J.EC.l)  DVX(  I  ,  J  )  - ( R A  *T (  I , J ) *S I N ( AF ) — 4 ) * HR  2-2  *  ZB (  I  ) * 

*  V  VX  ( I  *  J+  1  ) 

IF( J.EG.M  )DVX(  l  , J)  =  (RA*T(  I  , J ) *SI N { AF ) -4 > *HR2-2*ZB ( I  )* 

*  VVX(I,J— 1) 

q v  X ( I  ,  J ) = ( D V  X (  I  ,J)— AVX (  I  »J)*GVX  (  I— 1  ,J)  )  /  (  B  V  T  S  (  I  , J  )  AVX 

*  (  I  , J ) *PV  X( 

C I -1  ,  J  )  ) 

10  CONTINUE 

VVX ( M  »  J  ) =C VX ( M  ,  J  ) 

DO  4  11  =  2,  MD 

4  VVX (  I  ,J  )  =  GVX(  I  ,  J  3-FVX (  I  ,  J  )  *VVX (  I H 1  , J ) 

9  CONTINUE 
DO  a  J=  1  »  N I 

co  a  i=  l ,  m 


« 
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CI-0NEG1*(VVX(  I  ,  J  )  —  V  X  (  I  ,J)  ) 
VX {  I  . J )  —  V  X  (  I  , J ) +0  I 
a  CONTINUE 
RETURN 
ENC 


91 

C  SUBRCLT  I  N  E  SUET  IS  TC  SCLVE  THE  ENERGY  EQUATION  (2-10) 

C 

c 

c 

SUBROUTINE  SUET ( IT  *AT  ,  C  T  , PT , F  ) 

DIMENSION  VX(21,31),T(31»3l)»VR(31»3l)»VC(31»31)»8VTS( 

*  31.31). 

1 A  T (  3  1  .  3  1  )  .  C  T ( 3  1  ,31  )  ,  D  T ( 3  1  , 3  l  )  ,PT(31  ,3  1  )  ,GT(31  ,31  ),TT(3 

*  1.31). 

CCITA  (31  )  , R { 3  1  ),ZA(3l),ZE(31)  , ZC ( 3  1  ) 

INTEGER  44  M ,N ,M  I  ,  M  »  N  D  »  M D ♦ V S  TD  ,  NO 

COMMON  HR  ,  HR  2  ,  FC  ,HC2  ,  AF  ,  RAC  •  RA  ,  PR  ,  STD  •  OK  ,  ZA  ,Z8.ZC.  R  ,C  I 

*  TA.EVTS. 

1 V  X , T .VR.VC.GMEG1 .ONES 
COM  .VC  N  M  ,  N  ,  M  I  ,  M  ,  NO  ,  MO  ,  MS  TO  ,  NO 

T  T (  1  ,  l)  =  0*25*(TT(2»l  )  ❖  (  1-HR* VR (  t  ,1  )*PR/2)+2*TT(2,N/2M 

*  >+TT(2iNI)* 

C(  l +FR*VR ( 1  , l  )*PR/2  )  +HR2*VX (  I  ,  1  )-F) 

DO  S3  J-2.N 

TT(1,J)=TT(1,1) 

T  T (  1  ,  NI  )  =  T  T (  1  , 1  ) 

D  T  (  2 , 1  )•=  ( -2  )  *ZE  (  2)  *TT  (  2 , 2  ) -VX  (  2 , 1  )  *HR2-TT  (  l  .  1  )  *AT  (2  • 

*  1  )  +F 

CT { 2 , J  >  =  ( -  1  )❖ ( ZB (2  )  +ZC (2 ) *PR*VC (2, J ) ) *TT<  2 , J-l  )— ( ZB( 2) 

*  —  Z  C  (  2  )  * 

C  P  R  ❖  V  C  (2, J  )  )  m ( 2  ,  J  +  1  )-VX ( 2, J ) *HR2-TT(  1  .  J ) *AT( 2  ,  J )+F 
DT{2,NT)=<-2)*Z£(2)*TT<2,N)-VX(2,NI )*HR2-TT ( 1 , 1 ) ❖AT ( 2, 

*  N  I  )  +F 

Q  T  {  2  »  1  )  —  C  T  {  2  , 1  )ZBVTS(2,l  ) 

QT(2,J)=DT(2,J)/EVTS(2,J) 

G  T  (  2  ,  N  I  )=CT(  2»M  )/BVTS(2,NI  ) 

CG  99  1  =  3, M 

D  T (  I  ,  l  )={-2> *ZE ( I ) *TT (  I  * 2>-VX <  I ,  1  )*HR2+F 

DT (  I  . J  )=< ZC(  I  )  ❖PR ❖VC (  I , J  )  — ZD (  I  >  >  *TT(  I  , J+l )-<  ZQ( I  )+ZC ( I 

*  )*FP*VC(I,J) 

C ) *  T  T (  I  »J-1  ) - V  X (  I  ,  J ) ❖  H  R  2  4  F 
DT (  I  ,NI  )  -(-2  )*ZE(  I  )*TT<  I  ,N)-VX(  I  .NI  )*HR2HF 
G  T (  I  ,  1  ) - ( D  T (  I  .  1  )  - A  T (  I ,  1  )*QT( 1-1,1  )  )/{BVTS(  I  ,  1  )-AT 

*  (1,1  )* 

CRT {  I-  1  ,  1  )  ) 

QT<l,J)=(CT(I»J)-AT(I,J)*QT(I-l»J))/(DVTS(  I  , J  )  - A  T(  I  ,  J ) 

❖  ❖PI (  I -  1  ,  J )  > 

QT{  I  .  N  I  )  =  {CT(  I  ,M  )  —  AT  (  I  ,  N  I  )*GT(  I  —  1  .  N  I  )  )/(BVTS(  I.NI)— AT 

❖  (  I  ,  NI  )  * 

CPT (  I  -  1  ,  NI  )  ) 

99  CONTINUE 

T  T  (  M  ,  l  )  =  0  T  (  M  •  1  ) 

TT(M,J)-CT(M,J) 

TT  (  M  ,  NI  )  -  0  T  (M  ,M  ) 

DO  3  I  I  =  2  •  ND 
I  =  M  I-  I  I 

T  T (  I  ,  1  )  =  CT(  1,1  )-PT(I,l  >*TT(I+l*l  ) 


,  I 
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TT(I,J)=GT(I,J)-PT(I,J)*TT(I+1.J) 

3  T  T (  I  »  M  )  =  GT(  I  , N  I  }-PT(  I  ,N  I  )  *T  T  (  I  +  l  ,  N  I  ) 
£5  CC  NT  INUE 

DG  77  J-l  ,NI 
CO  77  I  =  l  *  M 

DT  T -C  NE  G  1  *  (TT (  I  ,  J  )-T {  I  ,  J  )  ) 

T (  L  ,  J  )-T (  I  »  J  )  40  T  T 
77  CONTINUE 
RETURN 
END 
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c 

c 

c 

c 

c 

c 


SUBROUTINE  SUE  VS  IS  TC  SCLVE  THE  VGRTICIITY  EQUATION 
ANC  MOMENTUM  ECUATICN  FOR  SECONDARY  FLC'rt  IN  EG.  (2-7) 
AND  (2-8)  USING  ECUNOARY  VORTICITY  M  ETHOO 


SUBROUTINE  SUEVS(S,VC,SS,AS»CS,PS*AVX,CVX,PVX) 
DIMENSION  VX(31.21),T(31i3l),VR(31,31)  • VC ( 3  1  •  3  1  )  « BVT S( 

*  31,31), 

2AS  <  3  1 ,3  1  )  ,  C  S ( 3  1  ,31  )  ,PS( 3  1  ,3  1  )  ,GS(3l  ,31)  ,  AVX ( 31  ,31  ),CVX 

*  (  3  1  ,3  1  )  , 

3FVX(31  ,31  ),QVC(31  ,31  ),V 0(31,31)  ,DS(31,31)  »SMI(3)*VMI(J 

^  )  , 

CCITA ( 3  1  )  ,R (3 1  )  . ZA (3  1  )  . ZE (3 1  )  , ZC<  3  l  )  ,DVG(  3  1 , 3 l ) , TRA (31, 

*  3  1) 

4  ,  S  S  (  3  1  ,3  1  )  ,S(  3  1 , 3  1  ) 

INTEGER* 4  M,N,MI  ,NI,ND, M C. MS TD, NO 

COMMON  HR,  HR2, HC,  hC2,AF,  RAC, FA, PR, STD, CK,ZA,Z8,ZC,R, Cl 

*  T A , EVT  S  ♦ 

1 VX ,T  ,VR  ,VC  ,CMEGl  .CMES 
COMMON  M  ,  N  ,  M  I  ,  M  ,ND,MD,MSTO,NC 
CO  1?  J  =  2,N 
CO  10  1=2,  M 

TRA  <  I , J  )=RAC*CGS ( AF ) *{ S  IN (Cl TA ( J)  ) *( T ( I  +  1  , J  >-T ( I -1  , J ) ) 

*  *{FR/2)+ 

CCOS (C  IT A  {  J  )  >  * ( 7 (  I  ,  J+l  )-T ( I , J-l  )  )*ZC (  I  )  ) 

10  CONTINUE 
K  =  0 

VM  I  (  1  )  =  0 .0 


20  K-K  +  l 

VO ( M  1  ,  J  )  =  VMI  ( K ) 

DO  <3  1  =  2  ,  M 

DVO(  I  , J  )=TRA (  I  , J  )-( ZB{  I  >  +ZC (  I  )  *VC (  I  , J )  )*VO(  I  ,  J-l  )-( ZB{ 

*  I  )  — 

CZC  (  I  ) * VC (  I  ,J))*VC(I,J  +  l) 

IF(  I  ,00,  M)  CVC(  I  ,  J )  =  T  R  A  (  I  ,J)— (  Z3(  I  )  TZC  (  I  )  *  V  C.  {  I  ,J)  )*V0( 

*  I  , J— 1  )  - ( Z8 (  I 

C)  —  ZC (  I  )  *  V  C  (  I  ,J)  )  *  V  G  (  I, J  +  l  )  —  C  V  X (  I  »J)*VC(MI  »J) 

9  CONTINUE 

GV0(2,J)  =  CV0(2,J  )/CVTS(2, J ) 

DO  8  1  =  3, M 

8  Q  VO  (  I  *  J  )•=  (  DVO  (  I,  J  )— AVX  (  I»J)*GVO(I-l»J)  )/(  EVTS(  I  ,  J  )  —  A  V  X 

*  (  I  ,  J  >  * 

CPVX (  I  -  1  ,  J  )  > 

VU(M  »  J)  =  CVC(M»J ) 

DO  4  11=2* MD 


* 


I  —  Ml  —  I  I 

V0(  I  ,  J)=QVO(  I,  J)-PVX(I,J)*VC(I  +  1,J> 

. . . . 

,„e,TS, 

*  p  s  <  i  - 1  •  J  >  > 


J ) 
J ) 
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2  QS(2»J)-CS(2»J)/EVTS(2*J) 

3  CONTINUE 

DS ( M I  *  J )={-l  )4(ZE(Ml)*(SS(MI»J-l)+SS(MI  ,J+1  )  )  )+HR2*VG( 

*  M  I  .  J  ) 

GS( M I  *  J )= ( DS ( M  I  , J  )-AS  <  M  l  , J ) *QS (M  *  J  )  )/( BVTSt  M  I  ,  J  ) -AS ( M I 

*  *  J ) *PS ( N  »  J  )  ) 

SMI(K)  =  CS<H.J) 

IF ( K-2 )  50*51 *52 

50  VMI(2)=-100.0 
GO  TO  2  C 

5  1  VMI (2  )  =  ICC*SM  T  (  1  )/{  SMI  (  2)-SM  1(1)) 

GO  TC  20 

52  SS( M  I  * J  )  =  SVI  ( K  ) 

CO  18  C  l-  1  ,  MO 

l-U  I-  I  I 

18  SS<  I, J)=CS(I  , J)-FS{  I  * J)*SS(  I  +  l  ,  J) 

17  CO NT  I  NUE 

DO  76  J  =  2  »  N 
DO  7  6  1=2  *  M  I 

CI=GNES4(SS{  I  ,  J  )  -  S  (  I  *  J  )  ) 

S {  I  » J  )  =  S(  1  »  J  )  +  D  I 
78  CONTINUE 
RETURN 
END 
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C  SUBROUTINE  SUEVRC  ID  TC  SOLVE  THE  SECONDARY  VELOCITY  U 

C  AND  V  FROM  THE  DEFINITION  OF  STREAM  FUNCTION 

C 

c 

c 

S U GROUT  INE  SUOVRCIS, ERROR) 

DIMENSION  VX ( 31  ,31)  ,T<  3  1  ,3  1  )  ,VR(31  ,31)  ,VC(3l.31),OVTS( 

*  3  1  ,31  )  , 

CCI TA { 31  )  ,R (31  ) , 2 A ( 3 1 ) , 2E (3 1 ) , ZC ( 3  l  )  , VVR ( 3  1  ♦  3  1  )  , VVC<  3  1  , 

*  3  1) 

2  ,  S  (  3  1  ,  3  1  > 

INTEC ER *4  M,N*Ml,NI,ND,NC,MSTD,NO 

COMMON  HR,  HR 2, HC, HC2, AF,F AC, RA, PR, STD, OK, ZA,Z6,ZC,R, Cl 

*  TA.GVTS, 

1  V  X  ,  T  ,  VR  ,  V  C  ,  C  N  E  G  l ,  C  ME  S 
COMMON  M  ,  N  ,  M  I  , M  , ND , MD , N S  TD , NO 
DO  10  I  =  2  ,  M 
R  H  6  =  E  (  I  )  *FC*6 
RH12=R  (  I  MFC*  12 

VVR (1,1  )  =  (E*S (  I  ,2  )-S (  I ,3 )  )/RF6 

VVR(  I  ,2)  —  (  — S(  I,2)+8*S{I  ,  3  )  — S (  I  ,4)  ) / RH l 2 

VVR (  I  ,NI  )-<S(  I  , N  D  )  — 8  *  S (  I , N )  )/RH6 

V V R  (  .1  ,  N  )  —  (  S  (  I  ,  N-2  )-8*S  <  I  ,ND  )  +S  (  I  ,N  )  )/RF12 

DO  1C  J=3,N0 

VVR  <  I  ,  J  )  =  ( S(  I  ,  J-2  )-a*S ( I  , J- 1  )  +  8*S(  I  ,  3  +  1  )-S(  l  ,  J+2 ) ) /RH  l 

*  2 

10  CONTINUE 
HR  L  2=HR * 12 
DO  <5  J  =  2  ♦  N 

VVC.  (  2  ,  J  )  =  (  10*S(2,3)-18*S(3,J)+6*S(4,J)-S(5,J)  )  /  HR  1  2 
VVC  (  N  ,  J  )  =  (  S  (  iV-3  ,  J  )-6*S  (M-2  ,  J  )  +  18*S  (  M—  l  ♦  J  )—  1  0*S  (  M  ,  J  )  —  3* 

*  S  (  N  +  1  ,  J  )  )  / 

C  H  R  1  2 

DO  9  1-3  »  M  C 

V  VC (  I  ,  J  )- ( S (  1+2  ,J  )  —  8*  S (  1  +  1  »J)+3*S(  I  —  1  ,J)— S(  l -2  »  J  )  ) /HR l 

*  2 

9  CONTINUE 

WR  (  1  ,  l  )  =  (€*S(2,N/2+l  )  — S  (  3  *  N/2+  1  )  )  / ( 6*HR ) 

DO  8  J- 1  ,NI 
VVC (  1  ,  J)  =  0  .0 
8  V VR (  1  , J ) -VVR (  1  ,  1  ) 

DO  7  1  =  2  ,  M 

VVC {  I  , 1  )=0 .0 
7  VVC (  I  ,N I  )  =  C  .  0 
S0  =  0 .0 
S  U  M  =  C  •  0 
CO  6  1=1, M 

DO  6  J=  1  ,  N  I 
DVR  =  VVR(  I  ,  J ) —  V  R (  I  *3) 

CVC  =  VVC(  I  ,  J  )-VC (  1,3) 

SD=SC  +  AES  (  DVRMAES(CVC) 

S  UM  =  S  UM  +  AS  S  (  V  VC  (  I  »  J)  )  +  A  0  S  (  V  V  !<  (  I  ,3)  ) 


*  j||* 


VK(I.J)=VVP(ItJ) 
VC (  I  ,  J )-V VC (  I  ,  J  ) 
6  CONTINUE 

EPROR=SC/SUM 

RETURN 

END 


. 
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C  SUBROUTINE  TANK  IS  TC  CALCULATE  THE  NUSSELT  NUMBER. 

C  FR  ICT  ICN  FACTCR 

c 
c 
c 

c  eukt: 

C  FRE  1  : 

c 

C  FRE  2  I 

C 

c  nu  l  : 

C  NU2  : 

c  sgt  ; 

c  sgvx: 

c  sip: 

c  svxp: 

c 

SUB P CUT  INE  T ANK { SVXP . SVXT ♦ SGT .SGVX ,BUKT  » F  RE  1  .FRE 2. UNI . 

*  UN 2.  PER A. 

CRFE 1  .RFE2  .RNU 1  .RNU2  » ARNU .STP .FRE3 ♦ AVEN , AVEF .RFE 12. RFE  3 

*  .  F  M  > 

DI MENS ICN  VX ( 2  1  .3 1  )  ,  T{ 3  1  .3  I  )  .  VR  (3  1 . 2  1  )  .  VC { 3  1  .  3 1  )  . BV T S ( 

*  31.31). 

2GT(31).0VX(3l).VXP{31).VXT(31).TS(31). 

CCITAI31  )  »P(31  )»ZA(31)  . Z B ( 3 1  )  * Z C ( 3  1  ) 

INTEGER**  M  ,N  .MI  , M  .ND.MC.MSTD.NO 

COM5MCN  HR.HR2.HC  .HC2.AF  .RAC.RA.PR.STD.CK.ZA.ZR.ZC.R.CI 

*  TA.0VTS. 

1VX.T.VR  ,  V  C  .  C  M  E  G  1  .  C  N  E  S 
COM MC  N  M  »  N  ♦  M I  . N I  »  NO . MD »  MS  TO . NO 
P 1=3 . 14 15927 
DG  1C  J- 1  ,N  I 

G  T  {  J  ) =— ( T { M— 3 . J >/4-4*T(M-2,  J)/3  +  3*T(MD,  J)-4*T(M!. J)  ) /HR 
GVX(  J  )  =  —  (  V  X  {  M  -  3  *  J  )  /  4  —  4  *  V  X  (  M—  2.J  )/3  +  3*VX(M!D»  J  )  —  4*VX(M»  J 

*  ) ) / FR 
10  CONTINUE 

SGT-G  T(1)+GT(NI) 

SGVX=GVX( 1 )+GVX(NI) 

DU  9  J-2.N.2 

SGT^SGT  +C-  T  (  J  )  *4 
5GVX=SC- VX  +GVX  (  J  )  *4 
9  C.CNT  1NUE 

DO  e  J=  3  , NC . 2 
SGT=SC-T  TGT  (  J  )  *2 
8  SGVX=SG VX  ICVX ( J ) *2 
SGT=PC*SGT/(3#FI ) 

SGVX=FC*SGVX/(3*PI ) 

DC  7  1  =  2. M 

V  X  T  (  I  )  —  VX(  1.1  )*T(I»1)+VX(  I.NI  )  *  r  (  I  iM  ) 

TS {  I )=T (  I  . 1  )  +T {  I.NI) 

VXF (  1  )  = VX  (  I  .  1  )  +  VX (  I  .  M  ) 

DO  t  J  =  2  »  N . 2 


E'ULK  TEMPERATURE 

FRCCLCT  ION  CF  FRICTICN  FACTOR  AND  REYNCLDS  NO. 
OF  1ST  EXPRESSION 

PRCCUCT ICN  CF  FRICTICN  FACTOR  AND  REYNOLDS  NO. 
OF  2ND  EXPRESSION 

NUSSELT  NUMGER  OF  1ST  EXPRESSION 
NUSSELT  NUMBER  CF  2ND  EXPRESSION 
GRADIENT  CF  TEMPERATURE  ALONG  WALL 
GRACIENT  CF  AXIAL  VELOCITY  ALONG  WALL 
AVERAGE  TEMPERATURE 
AVERAGE  AXIAL  VELOCITY 
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VXT(  I  )—  V  XT  (  I  ) +VX {  I  »  J  )  *  T  (  I , J ) *4 
TS(  l ) -T  S  {  I  )  +  T(  I  » J)^4 

6  VXP  {  I  )  =  V  XP {  I)+VX{  I  *  J  )  4  4 
00  5  J-3jNC»2 

VXT(  I  )  =  V  XT (I)+VX(I  *  J ) *  T (  I  *  J ) *2 
TS (  i  )-7 S (  I  )  +T (  l  ,  J  )  42 
5  VXP{ II=VXF ( I ) +VX< [ , J ) *2 
V  X  T (  l  )=HC*VXT (  I  ) / 3 
TS(I)=PC*TS<I)/3 
VXP (  I  )=hC* VXP (  I  ) / 3 

7  CONTINUE 
SVXT=  0 .0 
STP-C  .  C 
SVXP^O  •  0 

DQ  4  1=  2  *  N  |2 

SVXT=SVXT  +  VXT(.I)*R{  I  )  *4 
STP  =  S  TP  +TS (  l  ) *R (  I ) *4 

4  SVXP=SV/P+VXP(  I  )  $  R  (  I  )  *4  / 

DO  3  1  =  3. VD, 2 

SVXT = SVXT  *  V  X  T  (  I  )  *  R  (  I  )  *  2 

STP=STP+  T  S  (  I ) *R (  I  ) *2 

3  SVXP=S VXPTVXP ( I ) *R ( I ) *2 
SVXT=HR*SVXT*2/(3*PI) 
STP=hP*STP*2/(3*PI) 

SVXP=HR  *S VXP *2/ ( 3*P  I  ) 

BUK T=S VXT/SVXP 
FRE  1=  4*SC-V  X/SVXP 

FRE 2= 8* { l-SIN ( AF )*RA*STP/4) /SVXP 

F R E  3=  0  *  O/SVXP 

UN  1  =  2  *  S  V  X  P  *  SG  T  /  S  V  X  T 

UN2  =  S VXP* *2/SVXT-FP* (  1/BUKT) 

RERA=2*SVXP*RAC 

RFE 1  =  FRE  1/  16 

RFE2-FRE2/ 16 

RFE3=FRE3/ 16 

RNU  1-UNP1  l/4€ 

RNU 2= UN 2* 1 1/40 
ARNU  =  (RNU 1 +RNU2 )/2 
AVEN= (UN  1+  UN  2  )  / 2 
AVEF— { FRE 14FRE2  )/2 
RFE 1  2  =  ( RFE  1  +  RFE2 )/2 
RETURN 
END 


. 
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Table  11  -  continued  fRe  and  Nu  with  ReRa  =  4,000,  a  =  45°  and  F  -  0 

(c)  Pr  =  5  cont.  _ _ _ 

Ra  (fRe),  (fRe),  fRe/(fRe)n  (Nu),  (Nu),  Nu/(Nu) 
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